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ABSTRACT 


An autopilot for the U.S. Marine Corps’ ducted fan hovercraft is designed using 
optimal control theory. Single input controllers are designed to govern the vehicle’s 
roll rate and altitude rate. The gyroscopic coupling between the vehicle’s pitch and 
yaw dynamics is examined and a multi-input controller is designed. A computer 
program called OPTCON is developed to generate optimal feedback control gains by 
solving the discrete matrix Riccati equation. This program is for use on portable or 
home IBM compatible computers. Graphic plotting of the time-varying gains and of 


the system time response is available for both monitor and hardcopy output. 


THESIS DISCLAIMER 


The reader is cautioned that computer programs developed in this research may 
not have been exercised for all cases of interest. While every effort has been made, 
within the time available, to ensure that the programs are free of computational and 
logic errors, they cannot be considered validated. Any application of these programs 


without additional verification is at the risk of the user. 
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I. INTRODUCTION 


A. THE CONTROL SYSTEM DESIGN PROCESS 
Since the beginning of time, man has sought ways to control the laws of nature. 

From the simple float regulator developed by the Greeks in 300 B.C. [Ref. 1: p 3], to 
the amazingly complex space shuttle of the 1980’s, control systems span the range of 
mankind’s efforts to govern his surroundings. The challenge for a control systems 
engineer is to use his knowledge, skill, judgment, and experience to systematically 
develop a solution to any of a number of different types of control problems. There is 
seldom only one right answer to a control problem. In general, there may be several 
alternate solutions to the same problem and the final product will probably be a 
compromise between them. It remains the responsibility of the engineer to choose the 
“best” solution that meets the performance criteria specified by the user. So, how does 
the engineer know where to begin when he is given a set of performance criteria for a 
svstem ? There is no set procedure carved in stone. There are, however, a few broad 
guidelines that give the engineer a rough idea of the tasks which need to be 
accomplished in his quest to design an effective control system. These milestones are 
as follows : 

1. Define the system. 

2. Specify the desired performance of the system. 

3. Identify the constraints under which the system must operate. 

d 


Translate the information from milestones 1, 2, and 3 into a mathematical 
model that can be simulated on the computer. . 


5. Use the available tools to develop a control system which satisfies the 
performance specifications. 


6. Evaluate the control system design using computer simulations. 
7. Modify the design as required to better suit the application. 


8. Incorporate the control design in a prototype system to test the ability of the 
system to tolerate real world non-linearities and non-ideal conditions. 


9. Modify and optimize the design until a satisfactory control system is realized. 
The first milestone listed above is not always as simple as it appears to be. In 
fact, defining the limits of the system may possibly be the most difficult phase of the 


design. If the engineer does not expend considerable effort in the definition of the 


ie 


problem which he is to solve, he might end up solving the wrong problem. The 
engineer needs to include enough parameters in his system to accurately model the true 
system without becoming overburdened computationally. This is as much an art as it 
is a science. The engineer will probably need to make simplifying assumptions and 
approximations which tend to widen the gap between the performance of the model 
and the performance of the real world system. 

Defining the desired performance of the system may or may not be left to the 
discretion of the engineer. If a strict set of specifications is handed to him, then the 
engineer=has little choice but to satisfy those specifications or be able’to defend his 
claim that they can not be satisfied. On the other hand, there may be considerable 
leeway for the engineer to make sweeping changes in the control system and still satisfy 
the required specifications. Several tools are available to measure the performance ofa 
control system. The classical design engineer holds fast to such measures as the gain 
margin, phase margin, root locations in the S plane or Z plane, and bandwidth. The 
advent of optimal control techniques has placed emphasis on the minimization of some 
cost function as a means of measuring system performance. All of these techniques 
have their place in the realm of control system design and it is the mark of a successful 
designer that he can incorporate any or all of the tools when the situation dictates. 

The third milestone is an important yet often overlooked element of the design 
process. Constraints on the system may include any or all of the following : 

1. Monetary cost 
2. -Admissible control inputs 
a. Saturation limits 
b. Observability of the parameters required for control 
3. Physical limitations 
a. Size 
b. Weight 
c. Minimum and / or maximum velocities, accelerations, etc. 
d. Initial conditions 
e. Final conditions 
f. Sampling rate and processor speed for digitally controlled systems 

The mathematical model is the link between the real world system and the design 

tools which the engineer has at his disposal. In general, most physical systems which 


need to be modelled can be represented by some set of differential equations. The 
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physical laws of nature lend themselves nicely to approximation by linear ordinary 
differential equations with constant coefficients. Non-linearities and random processes 
are also quite prevalent in many systems and the effects of these phenomenon can 
greatly complicate the engineer's effort to model a system. That is why he must have 
the expertise and experience to know how to make assumptions and approximations 
which simplify the problem at hand to a point where he can use the available design 
tools. 

The next step is to use design tools to develop a control system which satisfies 
the desired performance specifications. It is at this point that two schools of thought 
begin to emerge. The classical school of thought focuses on such design tools as the 
Root Locus Plot, Nyquist Plot, Bode Diagram, and Function Minimization. The more 
daring school of thought centers its attention on the maximum principle of Pontryagin 
and the method of dynamic programming developed by Bellman. The advent of digital 
computers in the mid 1950’s made these more powerful design tools realizable since the 
amount of computation required by them was prohibitive if not impossible to do by 
hand. In any event, the design engineer has a myriad of tools from which to choose. 

Once the design of the controller is accomplished, the next step is to integrate the 
control system with the system model and then simulate the entire system to evaluate 
its performance. A very useful method to perform this evaluation is to study the time 
response of the output variables of the system. Such parameters as rise time, peak 
overshoot, and settling time are typical values to be noted. The digital computer once 
again ts a very useful means of obtaining such information rapidly. 

Even the best of control designers is not apt to hit a bullseye on his first shot. 
Control system design theory does not guarantee success on every try. The method of 
trial and error is one with which all control engineers are familiar. Modifications to 
the control system are inevitable. ; 

Once the controller design is proven in simulation studies, it is time to test it out 
On a prototype or small scale model of the actual system. This phase of design can 
become costly if the designer has not thoroughly tested his controller on the computer 
first. In this phase of design, the non-linearities and random effects which were 
ignored or approximated during the modelling phase become significant factors once 
again. Conditions which were assumed to be ideal in the model now become non-ideal. 
The evaluation process begins all over again as these new disturbances change the 


performance of the system. 
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The next step is obvious. After evaluation of the controller in a real world 
system, the need for further modifications is again probable if not inevitable. Changes 
must be made until a satisfactory controller has been designed which meets the desired 


specifications. 


B. AROD 

It is the goal of this thesis to complete the first seven steps of the nine step 
design process discussed in the previous section. The system chosen for this endeavor 
is an airborne remotely piloted vehicle (RPV) called AROD. The acronym stands for 
Airborne ‘Remotely Operated Device. The United States Marine Corps initiated work 
on AROD early in 1986 and is attempting to introduce the vehicle into the operational: 
Fleet Marine Force during fiscal 1987. AROD is a slow, low-flying ducted fan vehicle 
powered by a vertically mounted, two cycle, two cylinder gasoline engine which drives 
a three-bladed propeller. See Figure 1.1. 

The vehicle is 38 inches tall, 32 inches in diameter, and has a nominal weight of 
85 pounds. It presently has a payload capacity limited to a miniature television camera 
and a canister of fiber-optic cable. The Marine Corps plans to use AROD for short 
range reconnaissance and over-the-hill spy in the sky surveillance. The fiber-optic 
cable provides two-way communication between the vehicle and the ground based 
operator. The uplink communication will consist of control commands while the 
downlink will provide real time surveillance and on-board status information. 

The primary flight mode is low altitude hovering with the axis of the spinning 
propeller oriented perpendicular to the surface of the earth. In order to translate 
horizontally across the earth’s surface, the entire vehicle must be tilted towards the 
direction of intended movement. The mechanism by which AROD 1s tilted for such 
translation consists of four control vanes located in the airflow downstream of the 
propeller wash. One pair of control vanes is designated as the rudder. The other pair 
is designated as the elevator and is oriented such that its axis of rotation is 
perpendicular to the axis of rotation of the rudder pair. See Figure 1.2. All four fins 
assume dual responsibility for aerodynamic control in that they also serve as ailerons 
for AROD. The control vanes are actuated by model airplane servos. These servos 
are limited to a maximum deflection of = 30° and a maximum angle rate of 50°/sec. 
A maximum translational velocity of 30 knots ( 34 mph ) in a no wind condition is 
desired. The translational velocity of AROD is proportional to the tilt angle created 


by the rudder and elevator control vanes. 
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Figure 1.2 AROD Control Vanes 


Vertical flight control of AROD is accomplished through a throttle controller 
which incorporates the same type of model airplane servo used to actuate the 
aerodynamic control vanes. The throttle controller increases or decreases the engine 


RPM as required to raise or lower AROD vertically. 


C. THE PROBLEM 

AROD is interesting from the standpoint of a control system design for several 
reasons. Most significant is the phenomenon of cross-coupled dynamics between the 
pitch and vaw subsystems. In addition, AROD is a Multi-Input Multi-Output system. 


These topics are briefly discussed in the following sections. 
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1. Gyroscopic Coupling 
AROD’s propeller has a spin velocitv of 7200 RPM in the hovering condition. 
This creates a large angular momentum vector along the spin axis of the propeller. 
Thus, AROD can be thought of as a large gvroscope with its angular momentum 


vector oriented perpendicular to the surface of the earth. Consider a cylindrical rotor 


x — 
Le 
wy 
h 
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Figure 1.3 Spinning Rotor Orientation 


spinning about the x-axis with an angular spin velocity w. See Figure 1.3. Let the 
rotor be of mass, m, with moments of inertia I, , I , and I, about their respective 


axes. These moments are defined in the three equations below. 


I = fff dy? + 27) dv (1.1) 
memos 2-) dV (Ge) 
I = fff 6x? + y?) dv (1.3) 
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In these equations, 6 is the density of the rotor and dV is an incremental volume 
element of the rotor. In the case of AROD, it is assumed that the propeller is spinning 
with sufficient angular velocity that its dynamics can be approximated by the dynamics 
of a cylindrical disk having the same mass, m, radius, r, and thickness, h. The 


moments of inertia of the propeller then reduce to the following : 


ne 1.4 
ial: (1.4) 
SomGr + ee a 
= (1.5) 
m(3r7 + h?) 
= (1.6) 
; 12 7 


Notice that the moments of inertia about the y-axis and z-axis are equal to each other 
due to the symmetry of the problem. 
The angular velocity, @ , of the rotor induces an angular momentum vector, 


H,, defined by the following equation. 


sg (ea Oe (1.7) 
If a coupled pair of forces, F, directed parallel to the z-axis is applied to the the spin 
axis of the rotor at a distance, d, from the rotor’s center of gravity, as in Figure 1.4, 


then a torque, M , results. The torque is given by 
M=F xd (1.8) 


If the rotor were not spinning with angular velocity, @ , then this applied torque would 
result in rotation of the rotor about the y-axis. The angular momentum of the 
spinning rotor, however, results in quite a different response to the applied torque. 
According to Newton’s Second Law, an external force applied to the center of gravity 
of a rigid body results in a change in the velocity of that body. A corresponding 


change in the body’s momentum also results. The changes in velocity and momentum 





-— Figure 1.4 Spinning Rotor with a Force Couple Applied 


are in the direction of the applied force. This law extends into the realm of angular 
forces, or moments, and angular velocities. In short, an applied moment, M , results 


in a change in angular momentum. 
M =— = —— (cle?) 


Notice that the change in angular momentum is in the same direction as the applied 
moment. This is the key to gyroscopic precession. By vector addition of H, and M, 
it can be seen that a new angular momentum vextor, Hew» results. The new angular 
momentum is displaced by an angle, y, from the initial angular momentum vector, H,. 
This angular movement is called precession and it occurs at a precession rate, r, 


oriented as shown in Figure 1.5 and defined by Equation 1.10. 


iy 





Figure 1.5 Resultant Angular Momentum With an Applied Torque 


a dy 
: dt \ 


It can be shown that H,, M,, and r are always mutually perpendicular to each other 


[Ref. 2: p. 335], and that these vectors are related by the expression 
Vi =" ee H. (1.11) 


The handy mnemonic for this relationship is that “spin follows torque”. In this 


example, the spin vector, H_, , that results from the applied torque, MI , is rotated in 


new 
the direction of the torque. Notice that the magnitude of the precession velocity is 


directly proportional to the magnitude of the applied torque. 
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In the case of designing a control system for AROD, the preceding 
development is quite important. Recall that the propeller spins at a nominal velocity 
of 7200 RPM. The angular momentum of this high speed rotor has significant effect 
on the flight dynamics of AROD. In order to change the orientation of this angular 
momentum, as is required to accomplish translational flight, a considerable torque 
must be applied. This torque is produced by the four control vanes located in the 
propeller downwash. Note that the gyroscopic nature of AROD introduces cross- 
coupling of the pitch and yaw dynamics. For instance, when a pitching torque is 
commanded about the y-axis via the elevator vanes, AROD must first undergo an 
initial yawing motion about the Z-axis. Similarly, a yaw command from the rudder 
vanes results in an initial pitching motion about the y-axis. These cross-coupled 
dynamics must be considered by the engineer when designing the controller for the 
elevator and rudder vanes. 

2. Multiple Control Loops 

Another feature which makes AROD interesting for the control engineer is the 
Multi-Input Multi-Output (MIMO) nature of its dynamics. There are basically four 
subsystems which need to be controlled in order to make AROD fly. These 
subsystems are : 

eee Iwoulerate 

2. Rate of vertical climb 

3. Piteneangle 

4.~ Yaw angle 
Classical design tools such as Root Locus diagrams and Bode plots are not well suited 
for MIMO system design. Instead, the usefulness of these methods is primarily limited 
to Single-Input Single-Output (SISO) svstems. These systems are generally represented 
in terms of their S domain or Z domain transfer functions. The poles and zeros of these 
transfer functions determine how the time response of the system will behave. By using 
the graphical and analytic methods available through classical design theory, the 
engineer can generally place the poles and zeros of his SISO controller in such 
locations as to obtain an acceptable time response for the system. The complex 
interactions that typically accompany a MIMO system can become impossible to 
represent in terms of standard transfer functions. Thus, classical design methods may 
become powerless for some systems. By developing a state space representation of the 


system, however, the interactions can be accurately modelled. Optimal control theory 


Zk 


is founded on the state space representation of control systems and, therefore, it seems 


logical to pursue this theory for AROD. The basics of optimal control theory are 
presented in the following chapter. 
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If. OPTIMAL CONTROL THEORY 


A. FEEDBACK CONTROL 
1. Why Use Feedback ? 

Feedback control is familiar to engineers from all disciplines. In its simplist 
form, feedback control is nothing more than using the present condition, or “state”, of 
a system to influence its condition in the future. 

The advantages of state feedback control [Ref. 1: p.97] can be summarized in 
four points : 


1. Assuming that controllability and observability conditions are satisfied, the 
transient time response of the system can be easily controlled and adjusted. 


2. The sensitivity of the system to plant parameter variation is reduced. 
Rejection of disturbance and noise signals is improved. 


4. Steady state errors may be eliminated or reduced. 


These benefits are not free. The penalty for using feedback control may 
include disadvantages such as : 


1. System complexity increases because additional sensors mav be required to 
measure the feedback states. 


2. Sensors contridute to an increase in: 


a. Cost 
b. Size 
c. Weight 


d. Measurement noise 
3. Closed loop gain is generally lower than open loop gain. 
Despite these potential drawbacks, feedback systems are widely used in all engineering 
fields. 


2. System Classification 
Feedback provides a system with the ability to moniter and alter its 
performance. As the process advances in time, the system is apprised of the changes 


that occur in its states. This state information may be real time or may be delayed by 
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some finite interval of time. In general, systems may be separated into two categories 
according to the nature of the signals they process. These categonies are : 

1. Continuous time. 

2. Discrete time) 

An analog electrical circuit is one example of the first type of system. In this 
case, tne voltage and current signals assume values over a continuum of time. That is, 
given any two instants in time, the changing values of these signals may be distinctly 
measured. This is true regardless of how closely the two time instants occur. Such 
svstems are usually described by a series of differential equations. The Laplace 
transform is extremely useful in allowing frequency domain analysis and design of 
continuous time systems. 

A microprocessor based system 1s an example of the discrete time system. The 
clocked signals in this type of system are represented by a sequence of numbers. 
Typically, a sequence of sampled data results from measuring an analog signal at 
specific intervals in time. The time between measurements is referred to as the 
sampling interval, or At. The sampling frequency, f,, is simply the inverse of At. For 
an analog system with a Founer transform bandlimited to a maximum frequency, f 


max’ 


the Nyquist frequency, f,, is defined in Equation 2.1 [Ref. 3: p. 138]. 


i (2.1) 
A general rule of thumb for the design engineer [Ref. 4: p. 404] is to sample a system 
such that 


f = 10f (2.2) 


This guideline for selecting a sampling frequency is based on the following 
considerations : 


1. Most systems are not strictly bandlimited. The choice of a sampling frequency 
greater than the Nyquist frequency compensates for contamination by higher 
frequency disturbances [Ref. 5: p. 30]. 


bh 


The sampling frequency should be fast enough to avoid aliasing. This 
distortion is generated during the convolution reconstruction of a time signal 
from its Fourier transform [Ref. 3: pp. 135-137]. 
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Discrete time systems are represented by difference equations and the Z-transform is 
the mechanism by which these equations are analvzed. More details of the 
comparisons between continuous time svstems and discrete time systems will arise 
during subsequent discussion. 


3. System Structure with Feedback 
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— Figure 2.1 Basic Control System 





The basic form for any control system is illustrated in Figure 2.1. It is the 
responsibility of the design engineer to determine any or all of the items designated in 
this schematic. In this section, emphasis is placed on the feedback gain “black box” 
shown in Figure 2.1. The two methods most commonly used in practice to determine 
feedback gains are pole placement techniques and optimal control techniques. The 
element of trial and error is inherent in both methods. 

Pole placement techniques include analytical methods as well as the frequency 
domain methods previously mentioned. These methods are best suited to low-order, 
linear, time-invariant, SISO systems. Although these methods are extremely useful for 
certain problems, no detailed explanation of these classical techniques is included in 


this thesis. 


Optimal control theory provides an alternative to classical pole placement 
techniques. A primary advantage of optimal control methods is that feedback gains 
can be computed for a much broader range of control problems. Specifically, optimal 
control provides solutions for high order, non-linear, time varying, MIMO systems. 
Such systems are intractable with classical methods. In addition, optimal control 
affords the designer the option to specify a performance criteria which is not linked to 
such standard time domain criteria as rise time, percent overshoot, and settling time. 
For instance, using optimal control theory, the design engineer may compute feedback 
gains Which result in a system that responds in minimum time to a given command 
input. Selection of a different performance criteria might result in a system that 
responds with minimum energy expendiature, minimum fuel, or minimum deviation from 
the reference command. Sound engineering judgement and a thorough understanding 
of the system dynamics are prerequisites for effective application of optimal control 
theory. No guarantee is made that the feedback gains obtained by optimal control 
theory will result in acceptable system response. The designer should evaluate the 
system response and modify his performance criteria in order to achieve the desired 


output. 


B. SYSTEM DEFINITION 
1. Continuous Time Systems 


The foundation of a successful control system is an accurate model of the 
th 


plant which is being controlled. The state space representation of a general n~™ order 
continuous time system is described by the following matrix state equations 
x(t) = A(t)x(t) + B(t) u(t) (2.3) 
y(t) = C(t) x(t) + D(t) u(t) | (2.4) 
e(t) = x(t) - r(t) - (2.5) 
u(t) = F(t) {x(t) - r(t)} (2.6) 


where the definitions in Table 1 apply to a system with €control inputs 


and m measurable outputs. 
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TABLE | 
STATE SPACE DEFINITIONS FOR CONTINUOUS TIME SYSTEMS 


Term Dimension Definition 


x(t) (n X 1) State vector 
u(t) (€ x 1) Control input vector (0 < @ S n) 
y(t) (m X 1) Output vector (0<m Sn) 


r(t) (m X 1) Command input vector 

e(t) (m X 1) Error vector 

A(t) (n X n) Plant matrix 

B(t) (n X &) Control distribution matrix 
C(t) (m X n) Output distribution matrix 

D(t) (m x £) Feedforward control gain matrix 
F(t) (€ X m) State feedback gain yector 





In this thesis, a linear time invariant system will be assumed. This allows the 
time dependency of the process matrices, A(t) and B(t), and the measurement matrices, 
C(t) and D(t), to be eliminated. Because optimal control theory requires that all n 
states be available for feedback, the output distribution matrix, C, 1s set equal to the 
identity matrix, I. This indicates that m = n and that the state vector is completely 
observable. In addition, the feedforward control gain matrix, D, is assumed to be 


equal to the zero matrix. These assumptions lead to the following simplified state 


equations : 
x(t) = Ax(t) + Bu(t) (mn 
y(t) = x(t) (2.8) 
e(t) = x(t) - r(t) (2.9) 
u(t) = F(t) e(t) | (2.10) 
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Figure 2.2. Continuous Time System 


The realization of such a system is schematically illustrated in Figure 2.2. 
2. Discrete Time Systems 

Optimal control theory is applicable to the continuous time system presented 
in the preceeding section. The remainder of this thesis, however, will focus on the 
application of optimal control theory to sampled data systems. The motivation behind 
this effort is to develop an interactive, user-friendly software package that can be 
implemented on a microcomputer. The digital nature of sampled data systems make 
them ideal for analysis and design with these high speed computers. The theory for 
optimal. control of discrete time systems is well developed and closely follows the 
development for continuous time systems [Ref. 6]. 

As was noted earlier, many digital systems are the result of periodic sampling 
of analog systems. This fact makes it necessary to mathematically connect the two 


types of systems. For an analog signal that is sampled at the frequency, f., a discrete 
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signal value is measured every t = kAt seconds. In this notation, k is an integer time 
index in the range OS kS (N-1) where N represents the last sample period of interest. 
Letting the sample period be denoted as 


Ate oh OL 


and substituting a discrete approximation for the derivative in Equation 2.7, 


a SkT) ~ x((k+ 1)T) - x(kT) _ § (2.12) 
a 
yields the discrete state equation : 
X((K+1)T) ~ (1 + AT)X(kKT) + TBu(kT) (2a) 
The analytic solution for the discrete problem ts given by : 
x(k+1) = @M x(k) + FT u(k) (2.14) 
y(k) = x(k) (2.15) 
e(k) = x(k) - r(k) (2.16) 
u(k) =  F(k) e(k) (2.17) 
where ® and [I are defined as : 
® = Al (2.18) 
r= ff Atas (2.19) 


th order discrete time system 


The vectors and matrices which describe an n 
with € control inputs and m measured outputs are summarized in Table 2. A graphical 


realization of such a system is illustrated in Figure 2.3. 


TABLE 2 
STATE SPACE DEFINITIONS FOR DISCRETE TIME SYSTEMS 





" Term Dimension Definition 
x(k) (n X 1) State vector 
u(k) (€ x 1) Control input vector (0 < € Sn) 
y(k) (m X 1) Output vector (OFS im =n) 
r(k) (m X 1) Command input vector 
e(k) (m X 1) Error vector 
@M(k) (he) State transition matrix 
T(k) (n xX 8) Discrete Control distribution matrix 
C(k) (m X n) Output distribution matnx 
D(k) (m X £) Feedforward control gain matrix 
F(k) (€ X m) State feedback gain vector 


The computation of the discrete process matrices, M and I, from the 
continuous process matrices, A and B, is readily accomplished [Ref. 5: p. 37] on a 
digital computer as follows. 


Define an auxillary matrix, IT as 
II = ji eA at (2.20) 


AT2 A2T? Alger 1) 
+ toe t+ 


=IT+ — 
2! 3! (i+ 1)! 





The terms in this Taylor series expansion are computed until the result is within a 


specified degree of accuracy. It behooves the programmer to set a very small tolerance 
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CONTROLLER SYSTEM 





Figure 2.3 Time Invariant Discrete Time System 


on the difference between successive terms in the expansion since this calculation is the 
critical link between the A and B matrices of the continuous time system and 
the ® and [I matrices of the discrete time system. Note that this calculation need 
only be done once for a given system with a fixed sample interval, T. The link 1s 


completed by using Equations 2.21 and 2.22. 


@® = 1+ ATI (2.21) 


r = Y1B (2.22) 


The subroutine PHIDEL listed in Appendix C performs the calculations required to 
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solve Equations 2.20, 2.21, and 2.22. A tolerance of 10°" is used for the allowable 


error between the last term used from the Tavlor expansion and the first term not used. 


ou 


3. Constraints 
The system is now defined in terms of the ® and [ state space matrices. The 
next step in the design process is to identify any constraints under which the system 
will operate. These constraints are as unique to the problem at hand as 1s the system 
itself. Therefore, no detailed discussion on constraints is appropriate without first 


defining a specific problem. This is done in Chapter III. 


C. THE PERFORMANCE MEASURE 
1.- Quadratic Cost Function a 
The central theme in discrete optimal control theory is minimization of a cost 
function, J, defined in Equation 2.23. 


N-1 

J = e(N)He(N) + Y [e(k)Qe(k) + u'(k)Ru(k) | (2.23) 
k=0 

where 

1" Heme nerve oes 

e(N) = State vector at the terminal time 

e(k) = State vector at intermediate discrete times 

u(k) = Control vector at intermediate discrete times 

— H = Terminal state weighting matrix 

Q = State trajectory weighting matrix 

R = Scaler control weighting matrix 

N = Time index at terminal time 

t = Matrix transpose operator 


2. Regulators and Trackers 
[t is imperative to note here that the error state vector, e(k), in Equation 2.23 
may not be the same as the system state vector, x(k), that 1s presented in the previous 


section. The e state vector is usually formulated in one of two ways : 


1. If e(k) = x(k), then the problem is a ‘regulator’ problem. The objective is to 
drive the system states to the origin during the time interval (1,N). 


2. Ife(k) = x(k) - r(k), then the problem is a ‘tracking’ problem. The objective is 
to drive the system states to have minimum deviation from the command input, 
r(k), during the time interval (1,N). 
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In order to extend the regulator problem to the tracking problem, the command input 
vector, r(K), must contain the same eigenvalues and structure as the x(k) state vector. 
It is possible, in many problems, to structure an approximate r(k) vector so that the 
use of the regulator solution is allowed. In the case of the regulator, the control input 


signal is generated as shown in Equation 2.24. 


u(k) = F(k) x(k) (2.24) 


a : — 


In the case of the tracking problem, the control signal is generated from the error 


signal as shown tn Equation 2.25. 


u(k) = F(k) {x(k) - r(k)} (2.25) 


The comparison between these two types of systems is demonstrated in their block 


diagram structures as shown in Figure 2.4. 


3. Performance Weighting Factors 
The H, Q, and R weighting matrices are the parameters by which the design 
engineer shapes the solution of an optimal control problem to best suit the problem. 
There are no magic formulas which instruct the designer on how to choose the values 
of these parameters. Intuition, experience, and patience are the keys to successful 
design. It is in selecting values for these performance weighting factors that the 
process of trial and error enters the design process. There are, however, some 
restrictions and general guidelines that partially direct the efforts of a design engineer. 
First consider the restrictions. Both of the state weighting matrices, H and Q, 
must satisfy all of the criteria below [Ref. 7: p. 753]. 
1. Dimension is (n X n) 
2. Real 
3. Symmetric 
4. Positive semidefinite 
In addition, the designer should never allow both H and Q to be equal to the zero 
matrix at the same time. The resulting cost function would be 
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Figure 2.4 Comparison Between Regulator and Tracker System Structure 
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Nel 
J = Y u'(k)Ru(k) (2.26) 
k=0 


It is simple to see that the cost function will be minimized by setting u(k) equal to zero 
[Ref. 7: p 755]. This would be disasterous since there would be no command signal 
available to drive the system states toward the desired state. It is permissible to set 
either H or Q equal to the zero matrix provided that they are never both zero 
simultaneously. The (€ x £) control weighting matrix, R, must be positive definite in 
order to assure the existence of a finite control [Ref. 7: p. 754]. . j 

Although there is no requirement for H and Q to be diagonal matrices, the 
usual practice is to select non-negative values for their diagonal elements and to set all 
off-diagonal elements to zero. This technique eliminates all cross product terms in the 
cost function. For example, consider a second order system containing states x, and 
x5. If diagonal matrices are selected for H and Q, then only the ns and (zaps terms 
will be weighted in the cost function. There will be no consideration given to the x,x, 
Or XX, cross product terms. 

Assuming that neither H nor Q is the zero matrix, it is suggested that the 
elements of these weighting matrices be selected so as to normalize the values of the 
states which they multiply (Ref. 6: p. 32]. For instance, suppose that x, represents the 
RPM of AROD’s propeller and x, represents the angular displacement in degrees of 
the throttle servo. State x, 1s expected to have a nominal value of 7200 while state X5 
may have a nominal value of only 10. Assume that element q,, which multiplies (a 
and element q,, which multiplies (x5)? are both set equal to one in an attempt to 
weight the two states equally. In terms of the cost function, the RPM will be weighted 
much heavier (approximately 720° times heavier '') than the throttle servo position 
angle. An appropriate solution is to set q,, = ( 1/720)? if q.. = 1. Thus, each signal 
is given equal consideration in the cost function. 

Notice in Equation 2.23 that the terminal cost term is not included inside the 
Ds operator. This means that the H term is only counted one time and therefore can 
contribute only once to the overall cost. The state trajectory term and the control 
term, however, are counted N times. If no compensation is made for this disparity, the 
cost of an error in the terminal state is likely to not have any impact on the control 
solution. [t seems that an additional scaling is required on the weighting elements. If 


the summation in Equation 2.23 is to be made over 500 discrete time intervals, for 
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instance, the normalized elements of H should be multiplied by 500 in order to weight 


the terminal states on the same scale as the state trajectory and control effort terms. 


Alternately, the elements of Q and R could be divided by 500. It is the relative 


magnitudes, not the absolute magnitudes, of these weighting factors which shape the 


control solution. 


4. Specific Types of Problems 


Optimal control theory can be used to solve any of the following types of 


problems : 


1. Minimal time 


2. Minimal control effort 
a. Minimum fuel 


b. Minimum energy 


3. Minimal error from a reference 


a. Regulator 
b. Tracking 


c. Terminal state control 


Each of these problem types requires minimization of a unique cost function in order 


to generate the appropriate control solution [Ref. 6: pp. 30-34]. The cost functions 


which are to be explored in this thesis are listed in Table 3. 


Goal 


Regulator 
Tracker . 
Terminal Control 
Minimum Energy 
Minimum Fuel 


TABLE 3 


TYPICAL COST FUNCTIONS 


Cost Function Additional Explanation 


1 = ¥ x'(k)Qx(k) 
> = ¥ el(k)Qe(k) e(k) = x(k) - r(k) 


e(N) = x(N) - r(N) 


> u(k)Ru(k) + J, 


5 Y[u(k) | 


J 
J 
J e'(N)He(N) 
J 
J 


All ¥° are performed over the interval0 Sk S(N-1). 





567. 


Of course, the cost function in Equation 2.23 is a general form which embodies all of 
the specific cost functions (except for minimum fuel) contained in Table 3. Proper 
selection of H, Q, and R will allow calculation of feedback control gains which cause 


the system to meet the specified performance goal. 


D. CALCULATION OF OPTIMAL FEEDBACK GAINS 

The method of dynamic programming is the workhorse which permits the 
calculation of optimal feedback control gains. Developed in the late 1950’s by R.E. 
Bellman,. this design tool provides a closed form solution for the minimization_of the 
quadratic cost function for a discrete time linear system (Ref. 6: p. 84]. 

The procedure involved in calculating optimal control gains is unique in that the 
computation begins with the final, or N™, discrete time interval and progresses 
backwards in time to the previous interval of the system process. This procedure in 
‘negative time’ is possible because the calculation of the gains do not require any 
information about the state vector, x(k). The sequence of calculations continues in 
negative time until a separate gain matrix is determined for every discrete sample 
period in the time interval (0,N). The resulting time-dependent gain trajectory is 
usually stored in memory so that the control signal, u(k), may be computed. 

An interesting and very useful property of the gain trajectory is its tendency to 
approach a constant valued matrix, F.., under certain conditions [Ref. 4: p. 354]. This 
constant matrix is referred to as the steady state feedback gain matrix. The conditions 
necessary for F.. to exist are: 

1. The system is controllable. 

2. The ® andI Matrices are time invariant. 

3. The H terminal state weighting matrix is equal to the zero matrix. 

4. The Q trajectory state weighting matrix is constant. 
5. The R control weighting matrix is constant. 

6. The number of stages, N, of the process is large. 
It is possible for the feedback gain matrix to approach F., without satisfying the first 
three conditions above. When all conditions are satisfied, however, a steady state gain 
solution is guaranteed provided that N is large enough. Just how large N must be in 
order to allow the gain trajectory to reach steady state is determined by the slowest 
time constant in the solution of the discrete matrix Riccati equation [Ref. 5: p. 259]. In 
practice, trial and error is the the most expedient method to determine how large N 


must be in order to ensure a steady state gain matrix, F... 


aE 


A series of three recursive equations [Ref. 6: p. 83] is used to calculate the gain 
trajectory, F(k). It 1s convenient to introduce a negative time discrete index, K, which 


is defined as follows 
K=N-k (2.27) 


Since k varies from (0,N-1) as forward time progresses, K varies from (1,N) as negative 
time progresses. Equation 2.28 is the solution for the transpose of the optimal 
feedback gain vector at each discrete time step. This equation is the solution to the 
well known discrete matrix Riccati equation. Equations 2.29 and 2.30 are aunxillary 


equations required to complete the calculations. The recursive matrix equations are : 


F(K) = - (P' p(K-1) F + RY! Crt P(K-1) O} (2.28) 
P(K) = ® + T F(K) (2.29) 
P(K) = YK) P(K-1) ‘P(K) + Q + FYK)R F(K) (2.30) 


with ‘negative time’ initial conditions 


F'(0) = 0 (2.31) 
(0) = 0 (2.32) 
P(0) =H (2.33) 


While somewhat difficult at first glance, Equatians 2.28, 2.29, and 2.30 are ideal for 
iterative solution by a digital computer. These equations are solved in the main 
OPTCON program listed in Appendix B. The reader who is not familiar with 
OPTCON 1s encouraged to review the discussion of this program in Appendix A prior 


to continuing with Chapter III. 
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Il. CONTROL SYSTEM DESIGN FOR AROD 


A. OVERVIEW 

The purpose of this chapter is to use optimal control theory to design an 
automatic flight control system for the U.S. Marine Corps’ remotely piloted AROD. 
In their initial form, the equations of motion which describe AROD’s dynamic 
behavior are extremely nonlinear and present a formidable challenge to the control 
system designer. For this reason, the equations are first linearized about a steady state: 
hover condition. The restrictions and assumptions under which the linearized 
equations of motion are developed are summarized below: 

1. The mass of AROD 1s constant with time [Ref. 8: p. 11). 


2. The propeller angular velocity 1s constant. 

3. Perturbations from steady state hover are small. This restriction requires that 
AROD pitch, roll, and yaw angular displacements be limited to less 
than 15¢ (7/12 = 0.2618 radians) [Ref. 8: p. 41]. 

4. Steady state pitch and yaw rates are zero. 

5. Initial side velocities are Zero. 

6. Initial bank angles are Zero. 

7. Initial angular velocities are zero [Ref. 8: p. 45]. 


In-order to elucidate the equations of motion, Figure 3.1 is provided for 
reference. The axis system in Figure 3.1 is known as a body-fixed coordinate system. 
The body-fixed axis system is thought of as being rigidly attached to the vehicle so that 
any change in the orientation of AROD with respect to the earth-fixed axis system 
(X’, Y’, Z’) results in a corresponding change in the orientation of the body-fixed axis 
system (X, Y, Z) with respect to (X’, Y’, Z’). The angles », 8, and wy, commonly 
known as the Euler angles [Ref. 8: p. 25], are measures of the roll, pitch, and yaw 
angles respectively between the body-fixed (X, Y, Z) coordinate system and the earth- 
fixed (X’, Y’, Z’) coordinate system. The angular velocities p, q, and r in Figure 3.1 
correspond to the roll, pitch, and yaw rates respectively. 

The automatic flight control system is logically separated into three subsystems 


according to the simplified equations of motion. The three control subsystems which 
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Figure 3.1 AROD Body-Fixed Coordinate System 


are hereafter designed are : 

I. Roll rate controller. 

2. Altitude rate controller. 

3. Pitch angle and yaw angle controller. 
Each of these controllers is designed independently from the other subsystems. 
Therefore, any cross-coupling which may occur across the subsystem boundaries is not 
accounted for. Each of the following sections is devoted to the design of a controller 
for one specific AROD subsystem. A detailed design is first presented in Section B for 
the roll rate controller in order to demonstrate the iterative nature of the design 
process. Section C presents the results for the altitude rate controller. In the interest 
of brevity, only the initial trial run and the final solution for the altitude rate controller 
are presented. The coupled dynamics of AROD’s pitch and yaw is examined in greater 


detail in Section D. 
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B. AROD ROLL RATE CONTROLLER 
1. The Roll System 

The purpose for roll rate control of the AROD is to allow the remote pilot to 
command a desired rotation velocity about the vehicle’s longitudinal, or x, axis. Such 
movement allows the camera aboard the vehicle either to slowly scan a selected ground 
sector or to terminate the rolling motion so that a target of interest can be further 
studied. The nature of remote sensing requires that the vehicle respond rapidly to a 
roll rate command. When the remote pilot locates a ground target, he needs to be able 
to swiftly bring the vehicle to a zero roll rate condition with negligible overshoot. Such 
movement is commanded through a twistable handgrip control located on the pilot's 
console. It is assumed that this roll rate command, p,, is limited to a step input 
of 1 radian/second (57.3°/second). Although no time response criteria are specified by 
the Marine Corps, it is assumed for the purpose of this work that the following design 
specifications for roll rate are required : 


1. Zero steady state error for a step input is required. 
2. The two percent settling time, tyo,, is less than 1 second. 
3. No overshoot is allowed. 


The simplified equation of motion which describes AROD’s roll rate 


subsystem is given in Equation 3.1. 


"Ge 
il 
— 


5 (3.1) 


aia 


The aileron servo dynamics are modelled in Equation 3.2 as a second order plant with 
a natural frequency, @, of 2 Hz and a damping coefficient, C, of 0.707. 


6. = 2605, - w76. a wu, (3.2) 


a 


The definitions in Table 4 apply to Equations 3.1 and 3.2. In order to apply the theory 
of optimal control, a suitable state space representation of the roll rate system must 
first be developed. Figure 3.2 presents the state space signal flow graph selected for 
this subsystem. 


41- 


TABLE 4 


VARIABLE DEFINITIONS FOR AROD ROLL RATE iat pes OF 
MOTION 


Variable Definition Value Units 


Vehicle Roll Rate TBD radians/second 


Aileron Effectiveness -21.29 seconds” 


Coefficient 


Aileron Servo <= 300 radians 
Deflection Angle 


Aileron Servo. S 50°/sec _radians/second 
Deflection Velocity 


Aileron Servo _ | 0.707 unitless 
Damping Coefficient 


Aileron Servo Ze radians/second 
Natural Frequency 


Control Input TBD volts 
to Aileron Servo 





By designating the output of each integrator in Figure 3.2 to be a state, the following 
third order state equations are derived : 


Pp 
is. x= 6. (3.3) 
6. | 
. QO -21.29 0 0 
x= | 0 0 l Xe 0 u, (3.4) 
QO -157.91 -17.77 157.91 
u.= F{x-r} (3-5) 


a 


Assuming that a unit step roll rate is commanded, the command input vector becomes 


De l 
r= |6..] = ; (3.6) 
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where the subscipt c indicates a command input variable. Now that the roll rate 
system and its constraints have been identified, the next step is to use the OPTCON 
program to design a suitable roll rate controller. 
2. Roll Rate Controller Design 
a. Choosing a Sampling Frequency 

In order to determine an appropriate sampling rate for the roll rate system 
given in Equations 3.3 through 3.6, the bandwidth of the open loop system is first 
determined. The open loop transfer function for the roll rate system is given in 
Equation 3.7. - 4 


P(s) (157.91) (21.29) 
a (3.7) 
U,(s) s? + 17.775 + 157.91 


The open loop Bode diagram for this transfer function is shown in Figure 3.3. The 
negative 3 dB bandwidth of the magnitude curve is approximately 12 
radians per second or 2 Hz. Using the criterion discussed in Chapter 2, The Nyquist 
sampling rate, f is 4 Hz. In order to avoid aliasing effects and to ensure that the 
sampled system is a reasonable representation of the continuous time system, it 1s 
decided to employ a sampling frequency that is at least five times greater than f,. A 
comparison of the effect of using various sampling frequencies is given in Table 5. The 
cost function which is used to obtain the optimal gains for this comparison is included 
in Table 5 and is hereafter referred to as the baseline cost function. The column 
labelled “Number of Stages Required” in Table 5 refers to the number of discrete time 
stages that must elapse before the optimal gain vector reaches its steady state value, 
F... Notice that the magnitude of the steady state gain vector is related to the 
sampling frequency. In general, a faster sampling rate yields larger feedback gains. 
Also notice that the sampling rate affects the amount of real time that is required for 
the gains to reach steady state. For example, when f, is 20 Hz, it requires 1.60 seconds 
for F.. to be achieved. When f is 40 Hz, however, a total time of 2.23 seconds elapses 
before F.. is achieved. These considerations are important if the control gains are to 
be dynamically implemented. In the case of this design, the control implementation is 
limited to steady state values of the optimal gains. The unit step time response 
obtained by using steady state optimal feedback gains is observed to be nearly identical 


for all three of the sampling rates listed in Table 5. Specifically, the roll rate state, x,, 
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TABLE 5 


EFFECT OF SAMPEING FREGEENCY 
ON ROLL RATE SYSTEM 







Baseline Cost Function 











Sampling Number 
Run Ce Tr F of 
(Hz) 5 Stages _ | 
Required | 
BZ 
} 2 


time response exhibits an overshoot of approximately 3.7% and a 2% settling time of 
1.3 seconds. See Figure 3.4. This implies that any of the three sampling rates 
examined is acceptable. Because it is generally considered good practice to implement 
small feedback gains when possible, it is decided to employ a sampling frequency of 20 


Hz for the remainder of the roll rate controller design. 
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Figure 3.4 Roll Rate Time Response 
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b. Methodology 

At this point, four different groups of cost functions are examined. For a 
given cost function group, the H and Q matrices are held constant while the control 
weighting factor, R, 1s varied within the range (0.01, 100). Approximately 15 runs are 
made for each cost function group with a different value of R inserted into the cost 
function for each run. After the steady state gains are determined for a given run, they 
are implemented into the control equation and a time response for the roll rate state, 
Xy is obtained. The percent overshoot and 2% settling time are recorded for each Xy 
time response. A summary of this information is presented for the four cost function 
groups in Tables 6, 7, 8, and 9. Following each of these four tables, there appears a 
graph of two time response parameters, percent overshoot and settling time, versus the 

value of the control weighting factor, R. 
It is hardly necessary to include a time response graph for all of the 59 
total runs examined. It is instructive, however, to compare the time responses for a 
selected set of cost functions. Three runs in each of the parameter summary tables, 
Tables 6 through 9, are flagged with asterisks. These flags indicate that the roll rate 
time response graph for that run is included subsequent to the applicable summary 
table. Keep in mind that the criteria for the roll rate step response is specified to be 
such that there is no overshoot and the 2% settling time is less than one second. A 
discussion of the results from the four cost function groups follows the last figure in 


this series of tables and graphs. 


=— — 
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TABLE 6 
ROLL RATE PARAMETERS FOR GROUP 1 


Sampling Interval T = 0.05 seconds 


Run Control Steady State Percent Settling 
oe f See c Overshoot ene 
l Z 
0.01 -0.2944 . 4,03 Lote 
0.03 -0.2942 4.02 hea), 
0.05 -0.2940 4.02 ieS)l 
0.10 -0.2936 4.00 esi 
0.30 -0.2920 Broo es 
0.50 -0.2904 3.89 1.31 
1.00 -0.2866 3.74 1.30 
3.00 -0.2728 3.17. 1.30 
5.00 -0.2610 2.68 1.29 
7.00 -0.2509 P22 1.24 
10.00 -0.2380 1.65 0.87 
20.00 -0.2078 0.42 1.01 
30.00 -0.1886 0.00 ely 
50.00 -0.1646 0.00 1.47 
ee SLOORUO -0.1350 0.00 2.05 


* 


t 


I 
iz 
3 
4 
) 
6 
4 
8 
9 


See Figure 3.6. 
Seemrigiure 3.7. 
See Figtire acc. 
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TABLE 7 
ROLL RATE PARAMETERS FOR GROUP 2 


Sampling Interval T = 0.05 seconds 


Run Control Percent Settling 


Ce 


0.01 
0.03 
0.05 
0.10 
0.30 
0.50 
1.00 
3.00 
5.00 
10 7.00 
11 ** — 10.00 
12 30.00 
13 50.00 
14 *** 100.00 


Ds 


l 
z 
3 
4 
5 
6 
val 
8 
9 


See Figur 
Scesriour 
See Figur 


fy 


0.4803 
0.4786 


0.4769 - 
0.4728 


0.4576 
0.4441 
0.4159 
0.3442 
0.3024 
O20. 
0.2435 
Di le26 
0.1281 
0.0937 


e 3.10. 
e 3.11. 
e 3.12. 


Steady State 
G 


ains 
2 
-1.2269 
-1.2218 
-1.2168 
-1.2047 
-1.1598 
-1.1202 
-1.0381 
-0.8374 
-0.7251 
-0.6504 
-0.5735 
-0.3702 
-0.2969 
-0.2175 
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Overshoot 


4.37 
4.35 
4.33 
4.28 
4.08 
3.88 
3.49 
Ze 
ee 74 
0.44 
0.00 
0.00 
0.00 
0.00 


Time 
(sec) 


0.79 
0.79 
0.79 
0.79 
0.79 
0.79 
0.79 
0.73 
0.57 


_ 0.62 


0.70 
16 
1.46 
2.00 
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Figure 3.10 Roll Rate Time Response for Group 2 Run Number 4 
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Figure 3.12 Roll Rate Time Response for Group 2 Run Number 14 
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TABLE 8 
ROLL RATE PARAMETERS FOR GROUP 3 





Run Control Steady State Percent Settling 
in ei c Overshoot (es 
Z 
0.01 -2.6905 4.53 0.48 
0.03 -2.6141 4.54 0.49 
0.05 -2.5460 4.58 0.49 
0.10 -2.4028 4.65 0.49 
0.30 -2.0369 4.63 0.51 
0.50 - 1.8200 4.29 0.52 
1.00 . -1.5079 3.64 0.53 
1.50 -1.3280 2.87 0.54 
2.00 -1.2053 2.33 03 
3.00 - 1.0425 1.28 0.43 
5.00 -0.8576 0.00 0.50 
10.00 -0.6458 0.00 0.71 
15.00 -0.5427 0.00 0.87 
20.00 -0.4782 0.00 0.97 
30.00 -0.3986 0.00 1.16 
50.00 -0.3153 0.00 1.46 
17 *** 100.00 -0.2277 0.00 2.00 


* 


l 
2 
2 
4 
) 
6 
i 
8 
9 


See Figure 3.14. 
Sceceigure 3.15. 
See Figure 3.16. 





50 


(1) UOLOVA ONLLUDITM TOULNOD 
NT Oy or 





. oer e Fs ; 3 . 
OF: co a : 
pS a= ARS a4 ee ees 


ere ess ve 


Be ee RN in.) etei@inl a: oe) ete Wl cetmtamie ie neue 


al 
LOOHSH 


SNOLIONDA 1S00 & (NOD 


Sd LAW VaVd NOISHG ALVY 1104 








c 


E 


YALANVUVd ASNOdSAY AWIL 


v 


Figure 3.13 Group 3 Time Response Parameters 
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Figure 3.14. Roll Rate Time Response for Group 3 Run Number 4 
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Figure 3.15 Roll Rate Time Response for Group 3 Run Number 14 
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Figure 3.16 Roll Rate Time Response for Group 3 Run Number 17 
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TABLES 
ROLL RATE PARAMETERS FOR GROUP 4 


Sampling Interval T = 0.05 seconds 


Run Control Steady State Percent Settling 
Wecignt Gains Overshoot 


FY f 
0.01 -5.4007 9.57 
0.03 -4,3357 8.82 
0.05 -3.8586 8.67 
0.10 -3.2461 eG 
0.30 -2,3821 6.22 
0.50 -2.0316 5.96 
1.00 -1.6110 4.51 
3.00 -1,0739 1.33 
5.00 -0.8759 0.00 
10** 10.00 0.6551 0.00 
1 30.00 -0.4020 0.00 
12 50.00 -0.3175 0.00 
13 *** 100.00 0.2289 0.00 


* 


| 


l 
2 
3 
4 
5 
6 
a 
8 
3 


= | oee Higune dulic: 
Sue See Higtine oalo, 
whe Sete SANS 3740) 





64. 









(1) MOLOVA ONELUDIAM TOULNOD 
Ol Ol OT 








CVI 
LOOUSTEAAO | 
(ING: 


SNOLIONAA ISO) FP dno 


SUALAWVYVd NOISHG FTLVY TIOY 





0°8 09 Oe OZ 0°0 
YALINVAVd ASNOdSaY ANIL 


0°OT 


O°cl 


Figure 3.17 Group 4 Time Response Parameters 
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Figure 3.18 Roll Rate Time Response for Group 4 Run Number 4 
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Figure 3.19 Roll Rate Time Response for Group 4 Run Number 10 
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Figure 3.20 Roll Rate Time Response for Group 4 Run Number 13 
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c. Results 

(1) Group I Cost Functions. All three states are weighted with a value of 
unity for the Group I cost functions. This implies that the designer is attempting to 
minimize the error in all states with equal emphasis. Pursuit of this group of cost 
functions is made in order to determine general patterns of cause and effect. For 
example, it is apparent in Table 6 that the control weighting factor, R, significantly 
affects the magnitude of the steady state optimal feedback gain vector, F... By 
increasing the penalty on the control effort, the magnitudes of the feedback gains are 
reduced. Thus, if the control system is physically limited to some maximum value of 
control effort, then the Rterm is the logical parameter to adjust. The percent 
overshoot and settling time data from Table 6 indicates that R directly affects the time 
response as Well. From Figure 3.5, note that the value of R has negligible effect on the 
time response parameters for any value of R less than unity. Refer to Figure 3.6 in 
which R = 0.1. As the control weighting factor increases above unity, however, the 
tume response is dramatically affected. For values of R greater than 30, the time 
response exhibits no overshoot and the settling time appears to lengthen without 
bound as the system becomes increasingly slow. Refer to Figure 3.8 in which 
R = 100. Also notice in Figure 3.5 that there is no cost function in Group | that 
yields an acceptable time response which satisfies both of the roll rate critena. The 
cost function in Group 1 which yields a time response closest to the design 
specifications is run number 12 shown in Figure 3.7. This run is subsequently used as 
a basis for comparison of the time responses generated by the other three groups of 
cost functions. . 

(2) Group 2 Cost Functions. Because acceptable results are not obtained 
from the Group | cost functions, it is decided to place increased emphasis on the error 
in the x, state while reducing the emphasis on the error in the x, and x; states. This 
tactic is allowable because the maximum absolute values of the x, and x, states are 
significantly less than the constraints for the 6, and 5. servo states listed in Table 4. 
Table 7 summarizes the data for Group 2. Figure 3.9 evidences the relationship 
between the x, time response parameters and the control weighting factor, R. Notice 
in this figure that an acceptable time response is expected for any Group 2 cost 
function in which 10 S RS 20. Figure 3.11 shows the x, time response for run 
number 11 in which R = 10. This time response meets the required specifications for 


roll rate. Note, however, that the gains for this run are approximately 80% higher, on 
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the average, than the gains for the best run, number 12, in Group |. In order to 
reduce the gains so that only a small control effort is demanded, there is more work yet 
for bc come, 

(3) Group 3 Cost Functions. Making the penalty on the x, error 
state 10 times greater than the penalty on the x, and x, error states in the Group 2 
cost functions appears to be a reasonable mechanism for obtaining an adequate time 
response. In an effort to reduce the magnitude of the control effort, the ratios of 
yy B59, y,/M33, 4yy/Go9, and q,,/q3,, are increased to 100 for the Group 3 cost 
functions. Table 8 and Figure 3.13 present the data for this group. The time résponses 
obtained for this group are very similar to those obtained for Groups | and 2. Notice 
in Figure 3.13 that the overshoot is zero for all cases in which R 2 5. In addition, 
the settling time is less than one second when R S 20. The Steady state gains for run 
number 14 average only 42% greater than F.. for run number 12 in Group 1. Thus the 
hypothesis tested in the Group 3 cost functions 1s validated. 

(4) Group 4 Cost Functions. The cost functions in Group 4 penalize errors 
ony in the X, state and the control effort. That all other elements of H and Q are zero 
implies that no penalty is assessed against deviations in the x, and x, states. Table 9 
and Figure 3.17 present the data for this group. Run number 10 is deemed to be the 
most acceptable time response and is shown in Figure 3.19. While this design satisfies 
the design criteria, note that the steady state control gains average 93% greater than 
the most acceptable run in Group 1. This is the greatest increase in control gains that 
is observed. Also notice that the percent overshoot curve in Figure 3.17 increases 
upwards of 9%. This is much greater than the maximum overshoot of 4.65% observed 
in Groups 1, 2, and 3. For these two reasons, it is determined that the cost functions 
tested in Group 4 do not need to be further pursued. 

(5) Summary. Figures 3.21 and 3.22 summarize the information contained 
in Tables 6, 7, 8, and 9. It is interesting to note in Figure 3.21 that the first three’ 
groups of cost functions yield suprisingly similar curves for the percent overshoot of 
the roll rate system. That the Group 4 cost functions produce a much more erratic 
curve is attributed to the fact that no weight is placed on the x, or x, states in this 
group. The roll rate settling times in Figure 3.22 exhibit similar patterns for all four 
groups of cost functions. Note that in all cases there appears to be a minimum settling 
time possible when 2 S R S 10. For values of R > 10, the large emphasis on control 


effort produces small steady state gains which in turn yield a slow system. 
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d. The Final Design 
Based on the time response specifications and on the desire to design a 
control system which implements minimal steady state gains, it is decided that 
run number 14 in Group 3 is the best solution for a roll rate controller. The time 


response for this set of parameters appears in Figure 3.15. 


C. AROD ALTITUDE RATE CONTROLLER 
1. The Altitude Rate System 
“Because the primary flight mode for AROD is low altitude hover, it -is 
important that there be a reliable control system to maintain the vehicle’s vertical 
position relative to the earth’s surface. The throttle on AROD’s two cycle engine 
provides the mechanism for altitude rate control. Table 10 defines the terms which are 


involved in the altitude rate equations of motion. 


TABLE 10 


VARIABLE DEFINITIONS FOR AROD ALTITUDE RATE EQUATIONS 
OF MOTION 


Variable Definition Value Units 


Vehicle Altitude Rate TBD feet/second 


Engine Thrustto RPM _—0.0865 ft/seconds“/rad 
Dynamic Coefficient 


Change in TBD RPM 
Engine Speed 


Engine Lag 0.5 seconds 
Time Constant 


Engine Scale 837.8 rad/sec/rad 
actor 


Throttle Servo -S 300 radians 
Deflection Angle 


Throttle Servo | S 50°/sec _—radians/second 
Deflection Velocity 


Throttle Servo_ . 0.707 unitless 
Damping Coefficient 

Throttle Servo eo radians/second 
Natural Frequency 


Control Input TBD volts 
to Throttle Servo 
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By commanding a desired altitude rate, h >» the pilot sets in motion the following 
sequence of events : 

1. A throttle servo control signal, u,, is generated within the controller. 

2. The throttle servo position is adjusted. 

3. The actuator position commands a specific engine speed. 

4. A change in engine RPM causes a change in the vehicle altitude rate. 
The rate of change, h, of the vehicle’s altitude rate, h, 1s proportional to the change in 


engine RPM as shown in Equation 3.8. 


~~ 


fm . 


- @ 
i) 

= 

oO 


(3.8) 


where the dynamic constant, C,, is experimentally determined in wind tunnel tests. 


The engine is modelled as a first order lag system according to Equation 3.9. 


(-1/t,) 8, + (K,/t,) 5, (3.9) 


The throttle servo is modelled as a second order plant in Equation 3.10. 


§, = -26m6, - 075, + w7u, (3.10) 


The signal flow graph for this system is shown in Figure 3.23. The following state 


Space equations are used to design the controller for altitude rate : 


* 
h 
) 
x= : (3.11) 
9, 
5, 
0 0.0865 0 0 0 
e 0 -2 1675.5 0 0 
x = ct u, (3.12) 
0 0 0 l 0 
0 QO -157.91  -17.77 157.91 


74, 


SUINVNAG 
HOM 


SUOINVNAC 
ANISNA 





OAYAS 
ERP maOls| ape 


— Ee sees dei see eo See aS eS 


~. 


eyed 


a= 


Y3TIOULNOO 


Figure 3.23 Signal Flow Diagram for Altitude Rate Control 


) ae 


De ee (3.13) 


If a unit step is commanded for the altitude rate, then the command input vector 


becomes 


(3.14) 


2 2eo — 


2. Altitude Rate Controller Design 

The system given in Equations 3.12 and 3.14 is entered into the OPTCON 
program and a controller is designed according to a procedure similar to that explained 
for the roll rate controller. As before, a sampling rate of 20 Hz is used for all runs. 
Only two runs are hereafter presented because the lessons learned during the design of 
the roll rate controller apply equally as well to the design procedure for the altitude 
rate controller. The performance specifications for this control system are designated 
to be as follows : 

1. Zero steady state error for a step input is required. 
2. The two percent settling time, tyo,, is less than 5 seconds. 
3. No overshoot 1s allowed. 

__The first run is made using a baseline cost function. The results obtained for 
this run appear in Table 11. Notice that an incredibly large number of stages are 
required in order for F.. to be achieved using this cost function. If the gains were to be 
implemented dynamically at 20 Hz, more than 38 seconds would be required before the 
steady state gains are available. This clearly is not desirable since. the settling time 
must be less than five seconds. The unit step time response using steady state gains 
from this first run is shown in Figure 3.24. Even after 20 seconds, the desired altitude 
rate is not yet achieved. The cost function used to generate this solution is deemed to 
be unsatisfactory and a better solution is sought. 

The final run for the altitude rate controller is summarized in Table 12. The 
cost function for this run places 100 times more emphasis on errors in the x, state than 
on errors in the x, and x, states. The altitude rate time response shown in Figure 3.25 
exhibits acceptable performance characteristics. By choosing the cost function wisely, 


it becomes possible to design a satisfactory controller for this system. 
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TABLE 11 
INITIAL ALTITUDE RATE PARAMETERS 


Cost Function 


Steady State Number Percent Settling 

Gains of Overshoot Time 
Stages (sec) 
fy [5 f, fg Required 


-0.0544 -0.0461 -6.2776 -0.1948 763 0.00 > 20.0 


TABLE 12 
FINAL ALTITUDE RATE PARAMETERS 


Cost Function 


Steady State Number Percent Settling 
Gains of Overshoot Time 
Stages (sec) 
fy f, f, fy Required 


-0.3485 -0.0312 -4.5999 -0.1569 128 0.00 4.65 
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Initial Altitude Rate Time Response 
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Figure 3.25 Final Altitude Rate Time Response 
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D. AROD PITCH ANGLE AND YAW ANGLE CONTROLLER 
1. The Pitch and Yaw System 


Gyroscopic coupling between the pitch and yaw dynamics of AROD creates 


an interesting control problem. Refer to Table 13 for an explanation of the terms 


involved in the pitch and yaw equations which follow. 


TABLE 13 


-- VARIABLE DEFINITIONS FOR AROD PITCH AND YAW 
EQUATIONS OF MOTION 


Variable 


7G) ae 


c One ™& 


C4") 


| 


re 


4 O < 


Mme —& 
“ “ 


& 
~ 


wx 


Definition 


Vehicle Pitch Rate 
Vehicle Pitch Angle 


Pitch to Yaw 
Gyroscopic Coupling 


Elevator Effectiveness 
Coefficient 


Elevator Servo 
Deflection Angle 


Elevator Servo 
Deflection Velocity 


Control Input 
to Elevator Servo 


Vehicle Yaw Rate 
Vehicle Yaw Angle 


Yaw to Pitch 
Gyroscopic Coupling 
Rudder Effectiveness 
Coefficient 


Rudder Servo 
Deflection Angle 


Rudder, Servo 
Deflection Velocity 


Control Input 
to Rudder Servo 


Elevator’/Rudder Servo 
Damping Coefficient 


Elevator/Rudder Servo 
Natural Frequency 


80. 


Value 


TBD 
TBD 
-6.78 
-14.51 

<= 300° 

S 50°/sec 
TBD 
TBD 
TBD 

6.75 
-16.68 

S 300 

<= 50°/sec 
TBD 
0.707 
12.57 


Units 
radians/second 
radians 
seconds”! 
seconds”2 
radians 
radians/second 
volts 


radians/second 


radians 
-1 


2 


seconds 
seconds 
radians 
radians/second 
volts 

unitless 


radians/second 


The pitch and yaw equations of motion are given in Equations 3.15 through 3.18. 


q=Crt+M,6, (Sal) 
0=fqdt | (3.16) 

Eerie Oy (3.17) 

y = fr dt (3.18) 


Note that crosscoupling between the pitch and yaw equations enters via the two 
gyroscopic coupling terms, C : and C.. The values listed in Table 13 for these two 
coefficients are based on an assumed constant propeller velocity of 7200 RPM. 

As before, the elevator and rudder control vane servos are modelled as second 
order systems according to Equations 3.19 and 3.20. 


@o@ s 
= 2 2 
6, = -2Gwd, - w°6, + aru, (3.19) 
ee @ > 5 
6. = -20M6, - @°6. + w*u, (3.20) 


A coupled eighth order system results from Equations 3.15 through 3.20. The signal 
flow diagram which represents this MIMO system is given in Figure 3.26. 
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Figure 3.26 Signal Flow Diagram for Pitch and Yaw Angle Control 
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Defining the eight states to be: 


x ak q 6, 8 yw +r 6 5, (3.21) 


the state space equation for the pitch and yaw coupled system is defined as 


x = Ax + Bu (3.22) 
where 
0 -l 0 0 0 0 0 0 
0 0 -14.51 0 0 6.75 0 0 
0 0 0 0 0 0 0 
0 QO -157.91 -17.77 0 0 0 0 
A= (3.23) 
0 0 0 0 0 -!l 0 0 
0 -6.78 0 0 0 QO -16.68 0 
0 0 0 0 0 0 0 l 
0 0 0 0 0 O -157.91  -17.77 
0 0 
- 0 0 
0 0 
157.91 0 
B= (3.24) 
0 0 
0 0 
0 0 
QO 157.91 
and the multi-input control vector is 
=} ¢j = 3.25 
u = : eo Fea Th (3.25) 
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2. Pitch Angle and Yaw Angle Controller Design 
a. Methodology 
Notice in Equation 3.25 that the control input, u, isa(2 % 1) vector. Up 
to this point in the AROD control design, the control input has been limited to a 
scalar signal. The multi-input control that results from gyroscopic coupling between 
pitch and vaw requires special attention. Consider the following points :- 


1. The optimal feedback gain matrix, F, is determined from the solution of the 
discrete matrix Riccati equation. This solution requires that the inverse of the 
term (I! P(K-1) F + R} in Equation 2.28 be determined. - : ~ 


2. The cost function for a SISO system requires that the control weighting factor, 
R, be a scalar. 


3. The cost function for an n™ order MIMO system with @ control inputs requires 
that R be an (€ X @) matrix. 


Thus, for a SISO system, the solution for F is greatly simplified because the term in 
Equation 2.28 is a scalar quantity. For a MIMO system, however, a matrix inversion 
routine is required in order to solve for the optimal gains. Although computationally 
possible, it is decided for the purpose of this work that no matrix inversion routine is 
to be included in the current version of OPTCON. This implies that the ability of the 
OPTCON program to solve for optimal feedback gains is necessarily limited to SISO 
systems. This limitation is reasonable since a multitude of control problems can be 
reduced to single input systems. In the case of AROD, however, gyroscopic coupling 
iS a permanent feature of the pitch and yaw dynamics. Thus, a MIMO system is 
inevitable. The four step tactic used to aesee a control system for this non-trivial 
problem is as follows : 


1, First assume that the gyroscopic coupling terms, C, and C., are Zero so that the 
coupled eighth order system reduces to two independent fourth order systems. 


2. Use OPTCON to solve for the optimal feedback gains for the two independent 
systems. 


3. Implement the steady state gains so obtained for the fourth order uncoupled 
systems into a simulation model for the eighth order coupled system. 


4. Experiment with various combinations and modifications to the (2 X 8) 
feedback matrix, F SGT until a satisfactory time response is obtained for the 
pitch angle and yaw angle of the coupled system. 


Note that the design procedure listed above does not represent the most direct method 
to design a MIMO controller using optimal control theory. Rather, this method is an 


attempt to solve a complex problem using a tool that is designed to solve simpler 
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problems. For this reason, the results may not necessarily be expected to meet the 
same high standards required of the two previous control designs. The target 
performance specifications for the pitch and vaw control svstem are stated to be: 

1. Zero steady state error for a step input is required. 

2. The two percent settling time, tyo,, is less than 2 seconds. 

2 


Less than 10° overshoot is allowed. 


(1) Decoupling the System Equations. The decoupling procedure results in 


—_ 


two fourth order systems. The uncoupled pitch angle state space equations are: 


g 
q 

Xg = (3.26) 
9 
o; 

| 

0) -1 0 0 0 

a” 0 -14.51 0 0 

Xp = 0 0 0 Xg + ; u, (Ge27) 
0 OF 1579 =i] 157.91 

u.= F, {xg - rg} (3.28) 


If a unit step is commanded for the pitch amgle, then the pitch command input vector 


becomes 


a 90 
| 


ie (3.29) 


OnM™m CO DD 
oo CO = 


@ 
Q 
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The uncoupled yaw angle state space equations are : 


W 
ales 
r 
0 -| 0 0) 0 
_ e@ 0 0 -16.68 0 4 QO | 31) 
Sy. = x es 
nee 0 0 0 1 oe 0 c 
0 QO -157.91 -17.77 157.91 
Uae Xy - hy} (3.32) 


If a unit step is commanded for the yaw angle, then the yaw command input vector 
becomes 


oo € 
a 


Ty (3.33) 


onem 
-™” 
a 
oOo oO = 


The similarity between the dynamics of the uncoupled pitch and yaw systems is 
advantageous. For example, it is found that the elevator and rudder control gains, F, 
and F_., which are generated by OPTCON have identical steady state values. For this 
reason, only one set of steady state gains, F.., needs to be generated for each cost 
function. The individual elements of re are hereafter referred to as fs f, f,, and f,. 

(2) Solving for the Uncoupled Controller. The solution for F.. for the two 
fourth order systems is straightforward and follows the procedure established earlier in 
this chapter. Table 14 summarizes the data from the initial run. A unit step time 
response for the pitch or yaw angle controller implementing steady state gains from 
Table 14 is shown in Figure 3.27. Table 15 contains the data for the final run. The 
time response for this last controller appears in Figure 3.28. The steady state gains 
listed in Table 15 serve as the foundation upon which the coupled controller is 
subsequently designed. 
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TABLE 14 
INITIAL UNCOUPLED PITCH OR YAW PARAMETERS 


Cost Function 


Steady State Number Percent Settling 
Gains of Overshoot Time 
Stages (sec) 
fy ) f3 fy Required 


-0.1704 0.2417 = -0.2490 = -0.0858 96 : 4.15 


TABLE 15 
FINAL UNCOUPLED PITCH OR YAW PARAMETERS 


Cost Function 


=~ O&O O&O & 


Steadv State Number Percent Settling 
Gains of Overshoot Time 
Stages Sec) 
fy f, f, fy Required 


-0.3961 0.2808 -0.4158 -0.0259 38 : 2.50 
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Figure 3.28 Final Uncoupled Pitch or Yaw Angle Time Response 
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(3) The Coupled Feedback Matrix. In order to implement the SISO 
feedback gain vectors, F, and F_, into the coupled MIMO state equations, the (2 X 8) 


feedback matrix, F is formed as shown in Equation 3.34. 


mimo’ 


The O. elements of F_.__, represent those feedback elements which are not specifically 
generated by OPTCON. The success of the final control system is contingent upon 
proper selection of values for these elements of F_._.. For the purpose of the 
following discussion, the reader is encouraged to refer to the signal flow graph shown 
in Figure 3.26. : 

(4) Analysis. There are numerous ways to select the eight unspecified 
feedback gains in Equation 3.34. The seven schemes examined during the course of 
this design are summarized in Table 16. The first two columns in Table 16 identify the 
controller structure used to generate the elevator and rudder control signals, u, and u.. 
The third column refers the reader to the appropriate figure containing the pitch or 
yaw angle time response for that particular set of parameters. The last two columns 
summarize the time response data for each controller design. At the bottom of 
Table-16-are listed the exact numerical values for the controller gains. 

b. Results 

(1) Controller Number One. The feedback matrix, F nimo* in this first 
design assumes that the four states of the yaw equations have no influence on the pitch 
control input, ug. Similarly, the four states of the pitch equation have no influence on 
the yaw control input, Uy: As expected, due to the known coupling that exists 
between the pitch and yaw dynamics of the vehicle, the time response in Figure 3.29 
exhibits unsatisfactory performance. 

(2) Controller Number Two. For this design, the pitch angle state, x,, 


influences the yaw control input, u,,, while the yaw angle state, x,, contributes to the 


W’ 
pitch control input, ug. From the time response in Figure 3.30, it is apparent that this 


design is not satisfactory. 
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TABLE 16 
PITCH/YAW CONTROLLER DESIGN SCHEMES 


aa Time Percent Settling 
Design ae: Response Overshoot ‘Time 
Number Organization Figure (sec) 
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Figure 3.29 Pitch / Yaw Angle Time Response for Controller Number | 


92 


| (988) JW] IY 
00°OI 00°8 00°9 00°h 00°e 00°0 





On 0 00°0 


020 
AdQLISPUdL 4lels 


Ue | 


09°] 


ONTUS G4N1450 Y4SN 


Butzuawayduy , 


ASNOdSdd SWIL IX 


ee 


Figure 3.30 Pitch / Yaw Angle Time Response for Controller Number 2 
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Figure 3.31 


Pitch / Yaw Angle Time Response for Controller Number 3 
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Figure 3.32 Pitch / Yaw Angle Time Response for Controller Number 4 
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Figure 3.33 Pitch / Yaw Angle Time Response for Controller Number 5 


96 


REAL TINE (sec) 


(988) JWI TWIY 
00°h 00°¢ 00°e 00°] 00°0 





00°0 


Oh 0 


Oc'l 


09'T 


SNIUS G4NI440 Y4sn 
butyuawa]duy | 


ISNOdSIY WIL TY 


AdULISF tal dLluls 


Figure 3.34. Pitch / Yaw Angle Time Response for Controller Number 6 
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Figure 3.35 Pitch / Yaw Angle Time Response for Controller Number 7 


98 


(3) Controller Number Three. Because the gyroscopic coupling between 
the pitch and yaw equations is directly related to the pitch rate, x,, and the yaw 
rate, X¢, it is decided to include these two states in the makeup of the yaw control and 
pitch control respectively. This seems like a logical design tactic at first but the 
resulting time response in Figure 3.31 proves otherwise. This design is clearly unstable. 

(4) Controller Number Four. Notice in the coupled system signal flow 
diagram in Figure 3.26 that the coupling coefficients, C : and C., are nearly equal in 
magnitude but opposite in sign. This realization causes the designer to hypothesize that 


the sign of O1,, in F should be reversed from the value previously used in 


mimo 
controller design number 3. As shown in the time response of Figure 3.32, this 
technique yields promising results. Although a steady state error of 2.15% exists, there 
1s Merit in pursuing this design further. 

(5) Controller Number Five. By finetuning the value substituted into (1. 
in Fimo? the time response for pitch angle or yaw angle is made to satisfy the desired 
performance criteria. The time response in Figure 3.33 exhibits no steady state error, 
zero overshoot, and a settling time, t,o,, Of shghtly less than two seconds. This 
controller design, then, 1s completely satisfactory. 

(6) Controller Number Six. For this design, all eight states are allowed to 
influence both ug and Uy. The time response so obtained is shown in Figure 3.34. 
Even though the steady state angle is only 75% of the commanded value, this 
controller design appears to be potentially useful. By tuning the gains iteratively, it is 
believed that zero steady state error is achievable with this design. 

(7) Controller Number Seven. This final design is a modification to 
controller number 3. In this case, only the pitch rate and yaw rate contribute to the 
crosscoupled control signals. This design effort results in unsatisfactory performance 
as shown in Figure 3.35. 

c. The Final Design 

Controller number 5 is selected as the best design for a pitch/yaw angle 
controller. The time response for this design appears in Figure 3.33. Note that this 
design is based on feedback gains generated by OPTCON but that a modification to f, 
is required in order to obtain the final design. Thus, the controller is not optimal, by 
formal definition, even though optimal control theory provides the foundation for its 


development. 
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IV. CONCLUSIONS AND RECOMMENDATIONS 


A. CONCLUSIONS 
The lessons learned during the course of this work are as follows : - 


1. The cost function weighting parameters, H, Q, and R, play vital roles in 
determining the magnitude of the steady state optimal feedback gain matrix, 
F... These control gains, in turn, significantly affect the time response- of the 
controlled system. 


2. There is no magic formula to determine proper values for the weighting factors. 
A reasonable starting point is to use the baseline cost function in which all 
diagonal elements of H, Q, and R are assigned the value of unity and all off- 
diagonal elements are zero. 


3. The process of trial and error is prerequisite to the successful design of an 
optimal control system. Only through an iterative procedure does the engineer 
establish the true nature of the problem at hand. 


4. There are obvious trends to be aware of. These include: 


a. The sampling frequency, f., must be fast enough to avoid aliasing effects. 
As predicted, the use of a sampling frequency that is five to ten times faster 
than the Nyquist frequency seems to be adequate. 


b. As the selected sampling frequency increases, the optimal gains generated 
also tend to increase in magnitude. 
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The control weighting factor, R, for a SISO system can be used as a 
parameter to systematically alter the tume response of the system. As the 
relative weight on the control effort increases, the steady state gains tend to 
decrease in magnitude. This generally produces a slow system that exhibits 
little or no overshoot. On the other hand, if the value of R is decreased, 
the steady state gains can become so large that a very fast and highly 
oscillatory system results. : 


5. The controller design for a MIMO system is significantly more involved than 
the design for a SISO system. If the engineer can logically and accurately 
decouple the MIMO system into multiple SISO systems, then the design effort 
becomes much easier. As shown in the pitch and yaw controller for AROD, 
this method is feasible. 


B. RECOMMENDATIONS FOR FURTHER WORK 
The following areas present valid opportunities for useful expansion of this 


work : 


1. A parameter identification procedure which aids the design engineer in 
determining or estimating the A, B, ®, and I plant matrices is needed. The use 
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of a Fast Fourier Transform (FFT) algorithm is one possible solution to this 
requirement. Such a program could be used in conjunction with, but not 
necessarily integrated into, the existing OPTCON package. 


The OPTCON program is limited in that it does not generate optimal feedback 
gains for a MIMO system. A matrix inversion routine is needed so that the 
discrete matrix Riccati solution can be determined for any (n X ¢&) B or ® 
matrix. 


The theory of optimal control assumes availability of all states for feedback. 
The design process must account for the fact that all states are not always 
measured. In the case of AROD, the servo position and rate states are not 
aVailable for feedback. This means that an observer must be designed in order 
to provide the missing state information. 


The three control systems which are herein designed must be evaluated on the 
actual vehicle. Although computer simulations provide a wealth of insight, the 
proof of a good design rests in the ability of the system to function in the 
outside world. 
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APPENDIX A 
THE OPTCON PROGRAM 


1. OVERVIEW 

The purpose of this appendix is to describe in detail the OPTCON computer 
program which was developed in support of this thesis. OPTCON derives its name 
from OPTimal CONtrol. A previous edition of OPTCON by Professor H.A. Titus of 
the Naval Postgraduate School provided the starting point for the work that follows. 
The original OPTCON program allowed the user to input a state space system either in 
the continuous time domain or in the discrete time domain. Using matrix calculations 
to solve Equations 2.28, 2.29, and 2.30, this first version of OPTCOWN generated a table 
of feedback gains and immediately terminated execution. The motivation for 
improving the original OPTCON 1s fourfold. 


1. The design process is an iterative technique. The OPTCON program needs to 
be flexible enough to allow minor changes to be made to specific parameters 
without the requirement to re-initialize all of the cost function and system 
values. 


2. The gain trajectory table is not a convenient means by which to analyze the 
solution to an optimal control problem. A graph of the feedback gains verses 
time provides better insight. 


3. A time response of the state space is needed in order to allow the designer to 
_ quickly evaluate the performance of the system. 


4. The program should be user friendly. The original OPTCON demanded that 
the user flawlessly enter the correct response to every question on the first 
attempt. Woe be it to the user who accidentally types a letter in response to a 
question that requires a numerical answer. The program aborts and any effort 
that was spent in entering information is wasted. The frustration factor for 
such an unfriendly program is likely to leave the program sitting on the shelf 
with nobody to use It. 


With these points in mind, the OPTCON program is rewritten to provide an 
interactive, menu driven, user-friendly, optimal control design tool that capitalizes on 
the graphical capabilities of modern microcomputers. The program is written in 
MICROSOFT Fortran and is listed in Appendices B, C, and D. Appendix B contains 
the driver program, MAIN. Appendix C contains the majority of the subroutines 
which are called by MAIN. Appendix D contains the plotting subroutine, GRAPH, 
which makes use of the PLOT88 graphics software package. 
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In order to use OPTCON to its full potential, the user needs access to the 
following : 


1. A microcomputer capable of executing MICROSOFT Fortran based programs. 
During the development of OPTCON, an IBM AT computer configured with 
640 kbyte RAM and Intel’s 80287 math co-processor was used. 


2. Fortran, PLOTS88, and Math libraries. 

3. A monochrome or color graphics monitor. 

4. An Epson or LaserJet printer. 
Figure A.1 is provided to give the user a broad overview of the basic program -flow of 
OPTCON. The blocks outlined by solid lines represent program segments that must be 
performed during the initial execution of OPTCON. The blocks outlined by dashed 
lines represent program segments that are optional. The numbers that appear to the 
left of each block are referred to during in Section 2.d of this appendix. 

The remainder of this appendix illustrates the solution to a simple example 
problem using the OPTCON program. The intent here is not to focus on the specific 
example problem or on its solution but, rather, to focus on the capabilities and use of 
SPTCON. 


2, AN EXAMPLE PROBLEM 
a. The Second Order Integrator 
Consider the continuous time system shown in Figure A.2. The state space 
equations for this second order plant are derived by defining the output of each 
integrator to be a state. Using Equations 2.7 through 2.10 as a basis, the state 


e 0 1 0 | 
x(t) -|' HJ -[}| u(t) (A.1) 


are t A.2 
ea (A.2) 


equations become : 


w(t) = Fy(xp-ry) + f(xy - 19) (A.3) 
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Figure A.2 Second Order Integrator Signal Flow Diagram 


If the system is sampled every T seconds, Equation 2.20 yields : 


: “J 
ne (A.4) 
ot 


and Equations 2.21 and 2.22 yield 


el e A.5 

ane (A.5) 
T2/2 

r= [ra (A.6) 
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The preceeding calculations are done simply to allow verification of the PHIDEL 
subroutine in Appendix C. This subroutine converts an A and B continuous system to 
a @ and [ discrete system. For instance, assume that the system is sampled at 
f, = 100 Hertz. This means that T = 0.01. Equations A.5 and A.6 become 


p= |i 00 | a 

a | (A.7) 

ey a ie | | = a er 
01 


for the discrete time representation of the second order integrator. 
Before proceeding with the OPTCON program, the user is urged to verify that 
the system is controllable and observable. 
b. Controllability and Reachability 
According to Astrom, a system is controllable [Ref. 5: p. 104] only if “it is 
possible to find a control sequence such that the origin can be reached from any initial 
State in finite time.” Thus, controllability is a necessary condition for the regulator 
problem. A similar property called reachability is required for the tracking problem. A 
system is reachable [Ref. 5: p. 104] only if “it is possible to find a control sequence 
such that an arbitrary state can be reached from any initial state in finite time.” 
~—For continuous time systems, the properties of controllability and reachability 
coincide. That is, either a continuous time system is both controllable and reachable 
or it is both uncontrollable and unreachable. For a discrete time system, however, 
. controllability does not guarantee reachability. Reachability of a discrete time system 
does guarantee controllability. The reachability of a discrete system is important 
because the engineer should not spend a lot of time designing a controller that 1s 
physically impossible to implement. 
A simple test is performed to check the reachability of a discrete system. The 
(n X nf) reachability matrix, ie for an n™ order discrete time system with £ control 


inputs is defined as follows : 


w.=|T OF @T .. or | (A.9) 
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If the reachability matrix is of rank n, then the system is reachable. In the case of the 
example problem 


ee eae 
‘lor (A.10) 
Taking the determinant of W_ and setting the result equal to zero, the necessary 
condition for reachability is found to be that T # 0. Since an infinite sampling 
frequency is impossible to achieve, the system is reachable and it is reasonable to 
continue with the problem. . 

c. Observability 
In order to take full advantage of the optimal gain schedule, F(k), it is 

necessary that all of the states be observable. According to VanLandingham 
[Ref. 4: p. 308], a discrete time system is completely observable if it is possible to 
determine the initial state, x(0), of the system based on knowledge of the control input, 
u(k), and the output, y(k), over a finite number of time intervals. The test for 
observability closely follows the test for reachability. First, define the (mn *X n) 
observability matrix, W,, as 


e 
C® . 
w_ = | Co (A.11) 


cp(n-l) 


If the rank of Wis is n, then the system is observable. This implies that all of the 
states of the system are available for state feedback. If the system is not completely 
observable, then one or more of its states is not measureable. Either the system must 
Operate without the unobserved states in the feedback path, or an observer must be 
designed to estimate the unobserved states. [In the example problem, the observability 
matrix 1s 


(A.12) 


° 
oOo - Oo = 
_—_ ey —- © 
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This observability matrix is of rank 2 and the system is completely observable. Notice 
that if the output distribution matrix 1s 


c=[i 0] | (A.13) 


so that only x, is observed, the observability matrix becomes 


w =|) ° re 
a ol) T ee 


which is of rank 2 provided that T # 0. However, if only x, is observed, then 


Cc =[o | (A.15) 


oe aie A.16 
oT ou (A.16) 


and the system is not observable regardless of the sampling frequency. A state 
observer is needed to give an estimate of the x, state. 


d. Solution Using OPTCON 


_ This section is an introductory guide to OPTCON. The second order 
integrator problem is used to acquaint the user with the commands, features, and 
limitations of the program. The messages presented in this section are referred to as 
“screens” and are surrounded by numbered boxes. Neither the boxes nor the numbers 
by which they are referenced are actual features of the OPTCON program. They are 
simply used as devices to make the following discussion more understandable. 
Messages which are generated by OPTCON appear in standard print. Any responses 
which represent keyboard entry by the user are shown in italic print. If the response is 
to be y for “yes” or n for “no”, then either uppercase or lowercase letters are acceptable. 
If the response is to be an integer entry, as in the menu selections, the subprogram 
COMPARE is called to verify that the user has entered a valid integer. If the response 
is out of range of the acceptable values, or if the response is not an integer, then the 


program repeats the message until a valid response is entered by the user. 
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I. Starting the Program 
The user enters OPTCON by typing optcon on the command line. The 


following heading appears 


Screen 1 


OPTCON minimizes the following cost function: 
J =MIN ( X'(N) % H ® X(N) + Sum( X'(K) ®% Q * X(K) + U'(K) ® R *€ ULK))) 


The output of the program is the feedback gain matrix, F transpose, (F'), 
which, when multiplied by the State Vector (xX), 
yields a scalar control.(U). 


The following recursive equations were derived using dynamic programming, 
starting at the terminal time (N) and working backwards: 


(1) F°(K) ~(DEL**P(K~-1 J*PHI )/(DEL‘XP(K-1LIXDEL + R) F°(0)=0 
(2) PSI(K) = PHI + DEL¥F'(K) PSI(0 )=0 
(3) P(K) PSI'(K )*P(K-1LIMPSI(K) + Q + F'(K)RREF(K) P(0)=H 





2. Entering Initial Information 
The first entry required is a problem name. This name is used to identify 
the output file called OPTFILE which contains all matrices, gain trajectories, and time 


response trajectories for each run that the user requests during the problem session. 


Screen 2 


First enter the problem identification ( NOT to exceed 20 characters ). 


PROBLEM ID........second order example 





Next, the user selects the type of printer that 1s connected to the operating 
system. If graph hardcopies are not to be requested during the course of: the problem 
session, then the response to this question has no impact on the operation of 
OPTCON. If graph hardcopies are to be requested, however, then the answer to this 
question sets a flag that allows data to be properly formatted for the printer that is 
being used. Unpredictable results are expected if the user attempts to get graph 
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hardcopies from a printer that is not selected. In this example, the LaserJet printer is 
to be selected. 


Screen 3 


Select the type of printer that you are using 
( Answer lor 2 ) 


1} EPSON or THINKJET 
2} LASERJET 





Now the user enters the order, n, of the system. The maximum order 
which OPTCON can accept is eight due to the limitation of 64 kbytes per segment in 


the IBM AT microcomputer. The practice problem requires that a 2 be entered here. 


Screen 4 


Enter the ORDER of the system (up to 8). 2 





3. Entering the Cost Function 
Next, the number of discrete time stages, N, is entered. This number ts 
limited to 1000 due to the maximum dimension size of the arravs in OPTCON. The 


user should be aware of the relationship : 


tp = NAt (A.18) 
where 
t- = final time of the process 
N = number of stages over which the )) in Equation 2.23 
is to be performed. 
At = sampling interval 


In the example, let At = 0.01 seconds and t, = 10 seconds. This requires N = 1000. 
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Screen 5 


Enter the NUMBER of TIME INTERVALS (N) over which the cost function 
is to be minimized. (MUST NOT exceed 1000) 1000 





At this point, the weighting elements of the cost function are. to be entered. 
Assuming that the user wants to initially create a baseline solution for the problem, a 
reasonable starting point is to let all diagonal weighting factors assume a value of 


unity. The routine to enter the cost function begins with Screen 6. 


Screen 6 


Does cost function (J) include the State TRAJECTORY over all stages ? 
( Answer 1l»,2,0r 3 } 


1) 6YES...Set Q equal to the IDENTITY Matrix . 
2) YES...Each diagonal element of Q@ will be entered separately. 
3) NO....Set @Q@ equal to the ZERO Matrix. 





Selecting option | results in the Q matrix being echoed in Screen 7. The 
program then advances directly to Screen 11. 


Screen 7 


The states are weighted equally for the TRAJECTORY over all stages. 


The Q Matrix 


1.0000 .0000 
-0000 1.0000 





Selecting option 2 in response to Screen 6 allows the user to enter values 
for the diagonal elements of the Q matrix. All off-diagonal elements are automatically 
set equal to zero. For the sample problem, assume that the user wants q,, = 2.4 and 
qo. = 5. After entering the value for q,,, the user is prompted to enter the value for 
q55- Screen 8 results. 
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Screen 8 


Enter the elements of the Q matrix. 
(State weighting matrix for TRAJECTORY over all stages) 


Q(l,1) =? 2.4 
.a(2,2) =? § 





After the user enters all diagonal elements, the matrix is echoed in Screen 9. 
OPTCON then advances to Screen 11. 


Screen 9 


The Q Matrix 


2.4000 0000 
0000 5.0000 





Selecting option 3 in response to Screen 6 sets all elements of the Q matnx 


equal to zero. Screens 10 and 11 follow. 


| Screen 10 


The state TRAJECTORY is not included in your cost function. 


The @ Matrix 


-0000 -0000 
-0000 0000 
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Screen 11 involves a loop which allows the user to change any or all of the 
elements in the matrix that is currently being processed. This loop is subsequently 


referred to in this discussion as “the modify routine.” 


Screen [1 


Do you want to change any element of the matrix? 


-l)} YES...a SINGLE element. 
2) YES...the ENTIRE Matrix. 
3) 





Option | produces Screen 12 which allows the user to change a single 
element by identifying the row and column of the element to be changed. The row and 
column entries must be integers separated by a comma. Assume that the user wants to 


change q,. so that it equals 3. 


Screen 12 


Which element of the Matrix do you want to Change ? 
If I is the ROW and J is the COLUMN,....enter I,J 252 





The user is then prompted to enter the new value for the element that is to 


be changed. Screen 13 applies. 


Screen 13 
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If the user entered this situation directly from Screen 7 then the 


result 1s Screen 14. 


Screen 14 


The Q Matrix 


1.0000 -0000 
-0000 3.0000 


Any other changes? (Answer y or n) 





If the answer to Screen 14 1s y, then OPTCON returns to Screen 12 and 
allows changes to be made to individual elements of the matrix. Once the user is 
satisfied that that the Q matrix is correct, a 7 1s entered in response to Screen 14. At 
this point, OPTCON is ready to accept information relating to the H matrix. 


Screen 15 1s next. 


Screen 15 


Does cost function (J) include TERMINAL States ? 
( Answer 1,;2,or 3 ) 


1) YES...Set H equal to the IDENTITY Matrix . 
2) YES...Each diagonal element of H will be entered separately. 
3) NO....Set H equal to the ZERO Matrix . 





The program operation at this point is identical to the operation illustrated 
in Screens 6 through 14. The only difference now is that all of the matrix information 
applies to the H matrix. Assume that the user has set h,, = 5 and h,, = 15. 


Screen 16 results. 


Screen 16 


The H Matrix 


5.0000 -0000 
-0000 15.0000 


Any other changes? (Answer y or 
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Since this is the desired H matrix, a 7 is entered and the program advances 


to the section in which the user is asked to enter the value for R. Assume the desired 


value is to be 15.7. Screens 17 and 18 result. 


Screen 17 


Enter the value of the scalar R 
(Control input weighting factor) 


15.7 


Screen 18 


The scaler R = 15.7 


Any changes to R ? (Answer y or n) 





The program echoes the value entered and asks if there are any changes. If 


there are changes to be made, a y response returns the user to Screen 17. A nm response 
in Screen 18 indicates that the cost function, J, is now defined completely and Block 1 


of Figure A.1 is complete. The program advances to Block 2 of Figure A.1 
The user must now indicate if the problem to be solved is in the continuous 


e domain or in the discrete tume domain. Screen 19 applies 


Screen 19 


If you want to read in A and B matrices for a ‘CONTINUOUS TIME system 
@eeeaee#te@s#2#pefe@e#?t@8e@@eeees @ Ce ee ee Enter 0 


@#ee0eoente@#e?@e0¢08 @ 


If you want to enter PHI and DEL matrices for a DISCRETE TIME system, 
Matera tote cae © oiae see eee « Enter 1 





The sample problem is of the-continuous type and a 0 is the appropriate 


response to Screen 19. Screen 20 follows. 
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Screen 20 


You will enter the A and B matrices. 
Is this correct ? 





If a y response is entered in Screen 20, then the program advances to 


Screen 22. Otherwise, the message in Screen 21 appears. 


Screen 21 


You will enter the PHI and DEL matrices. 
Is this correct ? 





The program toggles between Screens 20 and 21 until the user enters y to 
one of these two options. Assuming that the continuous system is selected, the next 
section of the program allows the user to enter the A and B system matrices and the 


sampling interval, T. 


4. Entering the Continuous Time System Parameters 
The elements ~of the A matrix are sequentially entered as shown 


in Screen 22) 


Screen 22 


Enter the elements of the plant matrix--A. 
5 , 

i 

0 

0 


A(l1,1) 
A(1,2) 
A(2,1) 
A(2,2) 





Screen 23 echoes the A matrix and affords the user an opportunity to make 
any changes. The modify routine is entered unless the user responds to Screen 23 with 


a 3. In the case shown, all entries are correct and a 3 is appropriate. 
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Screen 23 


The A Matrix (Plant Matrix) 


.0000 1.0000 
-0000 -0000 


Do you want to change any element of the matrix? 
1) YES...a SINGLE element. 
2) YES...the ENTIRE Matrix. 





The elements of the B matrix are sequentially entered as shown 


in Screen 24. 


Screen 24 


Enter the elements of the distribution matrix--B. 


B(l,1) = 20 
Bikol) 2 7 1 





Screen 25 echoes the B matrix and once again allows the user to enter the 


modify routine if necessary. 


Screen 25 


The B Matrix (Distribution Matrix) 


-0000 
1.0000 


Do you want to change any element of the matrix? 
1) YES...a SINGLE element. 
2) YES...the ENTIRE Matrix. 
3) NO 


ANSWER. .......020008 


Since no changes are needed, the program now prompts the user to enter 


the sampling interval, T. The correct answer for the sample problem is entered 
in Screen 26. 


tl? 


Screen 26 


Enter the SAMPLE INTERVAL 





As usual, the response is echoed and the user is allowed to make changes 


until the correct value is entered. Screen 27 applies. 


Screen 27 


The SAMPLE INTERVAL OT = -0100 
Any changes to the SAMPLE INTERVAL ? (Answer y orn) 2n 





5. The Optimal Feedback Gains Calculated 

The program now has all of the information necessary to calculate the 
optimal gain schedule. The first step that OPTCON must perform is to convert 
the A and B matrices to the corresponding @andTI matrices. The subroutine 
PHIDEL in Appendix B performs this conversion. The resulting ® and I matrices are 
not displayed on the monitor. These two matrices are, however, recorded in the 
OPTFILE listing for the user’s convenience. If a discrete time system is initially 
selected in Screens 19 and 20, then the PHIDEL subroutine is bypassed. In either 
case, the gain schedule is now calculated using Equations 2.28, 2.29, and 2.30. This 
completes Block 3 of Figure A.1. As block 4 of Figure A.1 1s entered, the user may 


choose to view the gain schedule in tabular form. Screen 28 applies. 


Screen 28 


Do you want to see the gain schedule table on the screen ? 
(Answer yorn) y 





Since the user wishes to view the gain schedule table on the monitor, a y 1s 
entered in Screen 28. The user should remember two points before choosing to list the 
gain schedule on the screen: 


1. The gain schedule is automatically entered into OPTFILE regardless of the 
user’s response in Screen 28. If the user wants to record the exact values of the 
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gains, this output file may be examined later using the BROWSE, COPY, 
EDIT, PRINT, or TYPE commands in DOS. 


2. <A total of N lines of data are listed on the monitor when tabular output is 
requested in Screen 28. If N is on the order of several hundred or more, the 
design procedure is likely to lose momentum due to the lengthy delay involved 
in sending such a large array to the monitor. 


In order to illustrate the form of the data generated, Screen 29 lists a 
portion of the gain schedule table. Only the first ten time intervals are listed here for 
brevity. The actual sample problem lists a table with 1000 rows. 


Screen 29 


NEG REAL 
TIME TIME 
STEP INDEX F°(1) F°(23 


FOIE IEEE FEE IE EIEIO IE IEE TENE FIE FETE FE FETE FETE IEEE FETE IE HEPC FETE FE FE FETE FETE TEE 


SCWON ANS WD & 





Block 4 of Figure A.1l is now complete and Block 5 is initiated. The next 
option available to the user is to have OPTCOWN generate graphs of the optimal gain 
trajectories. If graphs are not desired, the user may answer 7 in response to Screen 30 
and the program advances to Screen 32. If plots of the gain trajectories are desired, 


then a y response is required in Screen 30. 


Screen 30. 


Do you want to see the gains plotted ? 


(Answer yorn) y 
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At this point, the program calls subroutine GRAPH. An internal flag is set 
so that the gain trajectory plots are sent to the monitor. A separate plot is generated 
for each gain trajectory. Thus, for an n™ order system, there are n separate gain plots 
produced. As each graph is generated on the screen, a pause is inserted so that the 
user may conveniently examine each one. Striking any key removes the current graph 


from the monitor and Screen 31 follows. 


Screen 31 


Do you want a hardcopy of this plot ? ( Answer yornm) y 





If a n is entered in Screen 31, then the program begins to generate the next 
gain trajectory plot for monitor output. By answering y in Screen 31, the user will 
automatically be provided with a hardcopy of the gain trajectory. A single hardcopy 
graph requires approximately 120 seconds on the Epson printer and approximately 90 
seconds on the Laserjet printer. Because of the superior quality of the graphs available 
from the later, all graphs contained in this thesis are generated on the Laserjet printer. 
As soon as the hardcopy is complete, OPTCON begins to generate the next gain 
trajectory plot for monitor output. It is important to note that the gain trajectories are 
plotted against the rea/ time discrete index, k. This means that the first gains calculated 
are those on the rightmost edge of the plot while the first gains implemented are those 
on the /eftmost edge of the plot. Thus, the term “steady state” as it applies to optimal 
feedback gains refers to the zero-slope property of the /eft side of the gain trajectory 
plot. The two gain trajectory plots for the example problem are shown in Figures A.3a 
and A.3b. When all n gain plots have been displayed on the monitor and/or have been 


printed as hardcopies, the program continues with Screen 32. 


Screen 32 


Do you want to change the NUMBER OF STAGES ? 


(Answer yorn) mn 





If the user is not satisfied with the initial choice of N, then a new value 
may be entered at this time by answering y in Screen 32. OPTCON presents Screen 5 


for this purpose and subsequently returns to the sequence beginning with Screen 28. 
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The most likely reason for the user to take advantage of the option in 
Screen 32 is that the gain trajectories do not reach steady state values in the allotted 
number of time intervals. By increasing N, the user may be able to force the gains to 
reach steady state. Since the gain trajectories in Figures A.3a and A.3b demonstrate 
steady state properties, there is no need to change N in Screen 32. 

6. The Time Response 

Block 6 in Figure A.1 involves calculation of the time response of the 
system based on the optimal gains computed in Block 3. The first option available to 
the user in this section is the phase plane graph of X, Verses X.. Screens 33 through 36 
represent the program sequence that results when a phase plane is requested with the 


following constraints : 


te = 10 seconds 

x,(0) = 0 = Initial condition on state x, 

x,(0) =O = __ Initial condition on state x, 

(1) =1 = Step forcing function on state x, 
(2) =0 = Ramp forcing function on state x, 


Screen 33 


Do you want to see a PHASE PLANE of X1 .vs. X2 ? 


(Answer yorn) y 


Screen 34 


For how many seconds ? 10 
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Screen 35 


the elements of the INITIAL STATE vector - XK(0) 


Screen 36 


Enter the elements of the COMMAND INPUT vector-R. 


R(1) 
R(2) 





The next option available is to select the method by which the optimal 


gains are to be implemented. Two choices are available. 


Screen 37 


Select a gain schedule....{ Answer lor 2 ) 


1) Use STEADY STATE gains over all steps . 
2) Use DYNAMIC gains . 





If the first option is chosen, then the state trajectories are calculated using 
Equations 2.14 through 2.17 such that the /ast gain matrix calculated, F(N), 1s 
substituted into Equation 2.17 during every cycle of the iteration process. The user 
must be aware of this procedure when selecting option 1 in Screen 37 because 
OPTCON makes no attempt to verify that the gains have indeed reached steady state. 
If the user selects option 1 when the gains do not exhibit steady state properties, then 
the solution is not optimal. If the gain trajectories do arrive at steady state prior to the 
N" stage, then selection of option 1 in Screen 37 may be appropriate. The time 
response obtained in this fashion represents the behavior of the system using a fixed 
gain feedback scheme. 
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The second option in Screen 37 causes the feedback gains to be 
dynamically implemented in the reverse order that they are calculated. The user ts 
cautioned that such implementation may not yield an acceptable time response. In the 
example problem, the gains reach steady state after approximately 500 stages. This 
corresponds to t, = 5 seconds when At = 0.01. Consider the case of a sampled system 
which has a transient time response longer than 5 seconds. The use of a'dynamic gain 
schedule would be disasterous in this situation. Because the gains progress towards 
zero as t approaches 5 seconds, the feedback channel is gradually eliminated from the 
system. The slow system, however, does not have enough time to reach steady state 
before the feedback gains go to zero. The error signal increases without bound and the 
system rapidly becomes unstable. Two simple solutions for such a situation are : 


1. Increase the number of time intervals, N. This causes the steady state portion 
of the dynamic gain schedule to become more predominate. 


2. Implement steady state gains instead of a dynamic gain schedule. 
After a gain schedule is selected in Screen 37, OPTCON begins to compute 
and save the state trajectories for x, and x, using Equations 2.14 through 2.17. The 


message in Screen 38 informs the user that the program 1s still executing. 


Screen 38 


Calculating Plotting Data 





After the state trajectories are computed, Screen 39 prompts the user. 


Screen 39 


READY TO DISPLAY DRAWING 
Strike any Key to continue. 





The monitor is cleared upon any keystroke and the x, verses x, phase plane 
subsequently appears on the screen. The graph remains on the screen until the user 
strikes any key. The monitor then clears and Screen 40 appears. 


> 


Screen 40 


Do you want a hardcopy of this plot ? ( Answer yorn) y 





If a nm is entered in Screen 40, then the program advances to Screen 41. 
Otherwise, the message in Screen 38 reappears. After a short delay, a hardcopy graph 
of the phase plane is automatically generated on the printer. See Figure A.4. The 


program then advances to Screen 41. 


Screen 41 


Do you want to see a time response of your system ? 


(Answer y orn) 





If a m is entered in Screen 41, then the first ran of OPTCON is complete. 
The program advances to Screen 44. Ifa y is entered in Screen 41, then the program 
prompts the user to enter parameters for the time response. Refer to Screens 34 
through 37. [t is not required that the user enter the same information for the time 
response that was entered for the phase plane. OPTCON recomputes the time 
response on every run. It is suggested, however, that the user carefully note the 
parameters that are entered for each run. Initial conditions and command inputs are 
recorded in the OPTFILE but this information does not appear on the graphs. After 
the gain schedule is selected in Screen 37, OPTCON begins to compute and save the 


state trajectories. When the calculations are complete, Screen 42 appears. 


Screen 42 


Do you want to see the time response table on the screen ? 


(Answer y or n) 
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Figure A.4_ Phase Plane for Second Order Integrator Example 
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A nm response causes the program to begin generating data for the time 
response plots. The user is cautioned that answering y in Screen 42 may result in a 
lengthy delay as the N rows of data are scrolled onto the monitor. The option to view 
this data on the monitor exists so that the user may gather exact numerical data 
without exiting OPTCON to examine the OPTFILE. A short segment of the tabular 
data appears in Screen 43. In the case of the example problem, this table continues 
until 1000 rows are displayed. 


Screen 43 


REAL 
TIME REAL 
INDEX TIME x(1) X(2) 


TIE ETE IE IEEE EIEIO EE EIEIO FE EHO IEEE TE TEE TE HEIE TE HE IE FE SE FETE FE NE FEE FETE IEE HEHE NE HEE TE TENE SEE FE EEE MEE ETE 


-0100 -0000 - 0000 
0200 .0000 .0099 
0300 .0002 ~0197 
0900 - 0004 -0292 
0500 .0008 0386 
0600 -0012 0479 
-0700 .0017 -0570 
0800 - 002% -0659 
.0900 .0031 0796 
- 1000 0038 0832 


1 
2 
3 
4 
5 
6 
7 
8 
9 
0 


po 





~— When the last row of the state trajectory table appears on the monitor, or 
if tabular output is not selected in Screen 42, then OPTCON begins to generate data 
for the state trajectory plots. Each state is plotted verses real time. During the 
calculations, the messages in Screens 38 and 39 prompt the user. State X, 1s plotted 
first and the n™ state is plotted last. The user may examine each individual graph on 
the monitor. By striking any key, the user clears the graph from the screen and the 
message in Screen 40 reappears. If the user does not desire a hardcopy, then a z 
response allows the program to process data for the next time response graph. Ifa 
hardcopy is desired, then a y is entered in Screen 40 and the procedure follows exactly 
as before. See typical time response plots in Figures A.5a and A.5Sb. After all n states 
are plotted, the program completes Block 8 in Figure A.1l. The first run of OPTCON is 
now complete and the user must answer y in Screen 44 in order to remain in the 
program. If an is entered in Screen 44, then execution terminates and the user 1s 


immediately returned to the DOS environment. 
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Screen 44 


This concludes the optimal control program (OPTCON). 


Do you want to rum the program again? (Answer y orn) y 





Assuming that the user desires to remain in OPTCON, a y 1s entered in Screen 44. The 


next section of the program is referred to as “the main menu” and is demonstrated 


— 


in Screen 45. 
Screen 45 


SELECT ONE OF THE FOLLOWING OPTIONS: 


1) Change the NUMBER of STAGES 

2) Change the TERMINAL state weighting matrix 

3) Change the TRAJECTORY state weighting matrix...Q 
4) Change the CONTROL weighting factor 


5) Change the present A and B matrices 

6) Change the SAMPLE INTERVAL 

7) Change the present PHI and DEL matrices 
8) Input an entirely NEW SYSTEM 

9) NO CHANGES...RUN 
10) EXIT the program 


SELECTION...( MUST be a number between 1 and 10 )...... 





It is mot necessary to describe in detail the operation of the main menu. 
The user should enter the integer value that applies to the particular modification to be 
made. If one of the first seven options is selected, the program responds as follows : 
1. Echo the current value(s) of the parameter(s) to be modified. 
2. Allow the user to keep or modify the selected parameter(s). 


3.° Return to the main menu for further modification, continued execution, or 
termination of the program. | 


If option 8 in the main menu 1s selected, then OPTCON returns to Screen 4 
and allows the user to enter new information for all parameters. In this case, no 
previous values are remembered by the program and execution proceeds just as if this 
is the first run. The OPTFILE, however, retains all information from any previous 


runs. 
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If option 9 is selected in the main menu, then the current values for all 
system and cost function parameters are written into OPTFILE to signal the start of a 
new run. Program execution begins with the gain calculation sequence represented by 
Block 3 in Figure A.1. Screen 28 applies. The user may rapidly skip through the 
intermediate steps of the program by answering v to several consecutive questions. 
For instance, suppose that the user changes a single parameter by selecting one of the 
first seven options in the main menu. In order to determine the effect of the changed 
parameter on the time response of the system, the following sequence of messages and 


— 


responses is used. | - 


Screen 46 
SELECT ONE OF THE FOLLOWING OPTIONS: 


1} Change the NUMBER of STAGES... ..cceccceseccseveNN 
2) Change TERMINAL state weighting matrix.....H 
3) Change TRAJECTORY state weighting matrix...@Q 
4) Change CONTROL weighting factor........++eeR 
5) Change present A and B matrices 

6) Change SAMPLE INTERVAL ask: ites tee ee coe el 
7) Change present PHI and DEL matrices 

8) Input an entirely NEW SYSTEM 

9) NO CHANGES...RUN 
10) EXIT the program 


SELECTION...( MUST be a number between 1 and 10 )......9 


Do you want to see the gain schedule table on the screen ? 


(Answer yorn)”n 


Do you want to see the gains plotted ? 


(Answer yorn) ” 


Do you want to see a PHASE PLANE of Xl .vs. X2 ? 


(Answer yorn) "nM 


Do you want to see a time response of your system ? 


(Answer yorn) yY 





At this point, the user may examine the system time response and evaluate 


the impact of the newly modified parameter. 
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By selecting option 10 in the main menu, the user is allowed to exit the 
program. The message in Screen 44 reappears as a safety mechanism to prevent 
inadvertant ejection from the program. If y is entered in Screen 44, then the main 
menu reappears and program execution continues. Otherwise, the program terminates 
and control is returned to DOS. 

7. OPTCON Summary 

The OPTCON program is designed specifically so that the user can easily 
modify problem parameters and rapidly obtain information about the effects of those 
changes. Tabular and graphical information is available both on the monitor and in 
hardcopy form. In an effort to make the program user-friendly, four techniques are 
employed : 

1. Menu driven options prevail. 

2. User input is screened for valid format. 

3. User inputs are echoed on the monitor. 

4. All data is written into an external file for later examination. 
The OPTCON program is quite useful as an interactive design tool for optimal control 
systems. Extensive use is made of its assets during the design of an optimal controller 
for the AROD in chapter three. 
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APPENDIX B 
OPTCON MAIN PROGRAM LISTING 


The following code is written in MICROSOFT Fortran and is intended to be 
used on an IBM compatible system. This is the main program for OPTCON and must 
be linked with the subroutines found in Appendices C and D. In addition, the Fortran, 
Math, and PLOT88 lbraries must be linked during the creation of an executable file. 


fj . 


NOfloatcalls 

NOdebug 

C 

C LL63.E0R LAST MOD 12JULY87 OK SDL 

c NEW selective state plotting 
: NEW state table formatting 


COMMON /BLK1/ A,B,PHI,DEL 
COMMON /BLK2/ BEGTIM, FINTIM,NPTS, 
+ XNAML , YNAML , PNAM1L , PNAM2L , PNAM3L 
COMMON /BLK3/ VTIME, VTIMSS , VY,VYSS, VXXSS , V&YSS 
COMMON /BLK4/ KFINAL,NSTAGE ,NSTP1,ORDERN, GNSKED , USERGN, FNEG, 
+ INPUT, Diy AVG, AVG2, MAXVAL , NINPTS 
COMMON /BLK5/ XNAME, YNAME , PNAME1, PNAMEZ2, PNAME3 
INTEGER*2 OPTION, ORDERN, IGOOD , CODE, NSTAGE, NSTP1,KFINAL,KPRIME, 
GNSKED, ‘NPTS, IOPORT, MODEL, ‘XNAML , YNAML , NCHAR1, NCHAR2, 
NCHAR3, STVAR, Ty, SKIP, OK, SYSTEM,GAIN, DTFLAG, PLTYPE, 
eaneae SCREEN, NINPTS , NINPP1, ORDNP1, GAINCH, GNSKD3, 
TPLO FLOTH 
REAL*4 psr(s, 8 ) 
FM(8) 
Er 8 


t++++ 


t+tttatt 
td 
a 
eH 
t 
10 


288) i DELING (8. /1), INPUT(8,1) 
D 8), FNEG(1000,85, 
TEMP , W(3 8). DT, TIME , PNAMIL, PNAM2L 
| apenas /VX¥SS(9) "AVG(BY. ,AVG2(8), 


oOo 
OOo 


CHARACTER*2 TEMP 

CHARACTER*3 ANS 

CHARACTER*20 NAME 7 
CHARACTER*30 XNAME , YNAME 
CHARACTER*51 PNAME1, PNAMEZ2 , PNAME3 


CHARACTER*5 OEM 
CHARACTER*4 HDG2(8) 
HDG(1 = trie) 
HDG(2 Se 
HDG(3) = 'F''(3)' 
HDG(4) = 'F''(4)' 
HDG(5 = URIS) 
HDG(6 = Pl 6) 
HDG (7 = fF os 
HDG(8 = 'F''(8)' 
HDG2Z(1) = 'X(1)' 
HDG2Z(2) = 'X(2)' 
HDGZ(3) =s2x(3)- 
HDG2Z(4) = 'K(4)' 
HDGZ(5)) = 
HDG2(6) = 'X(6)' 
HDG2Z(7) = 'X(7)' 
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HDG2(8) = 'X(8)' 
OPEN(9,FILE='OPTFILE' ,STATUS='NEW' ) 


© 
CRAKKAKAKAKKKAAAKAKKAKKKKAKKARKKKKKKKAKKKKAKKKKAAAKKKAKKKARAKARKKKKKAKKKAKKKKAKKK 


CHARAAAAAAAAAAAAA PRINT OPTCON HEADING and INPUT PROBLEM ID *xxxxxAAAKAA 
CAKKAAAKAKAKKKAAAKRKKKKKAARAKKKKKARKERKAKKKKKKAKRRKRKKKKKKARERKAKKKKARAARKAKAKKAKKKARA 


Ec 
Pe eae 
WRITE(*,2010 
PAUSE 
Bee oe 
READ (*,2020,END=1530 )NAME 


AKRAKKRKAKKAKKKRRAKKKKKKRAKAKKKKKRKARKKKKKEKKRERKAKKKKKAREAKRKKKKKAKEKAAKKAKKKARAARKK 


CHAAAAKAAAAAAAARA HEADING INFO FOR THE OUTPUT FILE KARRARR A ARAR 
CARAKKAKKAKKAKAKKKKKKKKRKKKKKKEKKERKERKEREKKRERKRKERKEKRKEKKERKKARERKRKRKRKRKRKRA RA 


WRITE(9, 2030 
WRITE (9,2040 
WRITE(9,2050 )NAME 
WRITE (9,2030 


CREKKRAKKAKAKKKKKRKRRAKKKKKKAAEAKAKKKRKKERAKKKKARAKRRKKKKKKRERKAKRKAKKKKRKKARAEKRKKARKAKKRKK 


CAAKKAKKARKRKARKARER ENTER PLOTTER/PRINTER MODEL TYPE AERKAARRARKK 
CRARAKKKKAAKKKAAKKAAAAKKRRRAKKRRAKKARRAAKKRRAKKKRRARARRRARKERARRAKAAKRRRKK 


5 WRITE(*,2055) 
READ (*,2070)TEMP 
CALL COMPARE (TEMP,1,2,CODE,IGOOD 
TE (CODE -£9-0)GOTO 5 


IF(IGOOD .EQ. 1)THEN 
IOPORT = 0 
MODEL = 1 
ELSE 
IOPORT = 0 
MODEL = 60 
ENDIF 
RARE ARERR RAEEERE RRR RRR ARERR RAKE RIOR RII CIR 
CRRA A A RAI INITIALIZE B ,DEL,USERGN MATRICES KKKKKAAKAKAR 
CRIA PKK KAKA KAKA RIK KAKA KR AKAE RRR KKAKRKRAKK AK KEKKARARKAKRAKKAARKRARR AK 
C 
DO 61 =1, 
p06 J =1,8 
B(I,J) = 0.0 
DEL(I,J) = 0.0 
USERGN(I,J) = 0.0 


6 CONTINUE 


c : 
CRAKKKKAAKKAKKARKAKKKAEKKKAAKKKKKKKRKKKKAKKKKKAKAKKKKREKKRAKAAKKKAKEKKKRAKKKKARAAKK 


CARKKAKKAKKKKKKAKK ENTER THE ORDER OF THE SYSTEM Kee KKK KKK 
CARKKAAKKKKKKKKKKKKKKKKAKKAKKREKARAAAAAKAKAKAKAKAKKAKAAKKAKKKAKKAKKKKAKKKKAKKKKKKAKK 


GS 
CRARKKAAKARKAKKK RESET FINAL ,GNSKED,GAINCH,GNSKD3 RARKKKAKKARK 
10 FINAL = 0 
GNSKED = 1 
GAINCH = 0 
GNSKD3 = 0 
Boro oor 
READ (*,2070) TEMP 


CALL COMPARE(TEMP,1,8,CODE,IGOOD) 
IF(CODE.EQ.0)GOTO 10 

ORDERN = IGOOD 

ORDNP1 = IGOOD + 1 


C 
CARKAKAKAKKKKAKKKRAKAKKKKKKKKKKKRKKKRKRKRKARRRRERERERKRRRRRKERKREREKRKRERKEREKRKKKKAKKK 


CAAAKARARARAKRARAR ENTER THE NUMBER OF CONTROL INPUTS RERRRKARAKKK 
CRAAAAKKAAAKKKKAAKAKKKARAKRRAAKKRERKRKRAKAKR AER KAKRARKKKRERKAKRRARAKAR RRA 
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15 WRITE (*, 2075) 
D (*,2070)TEMP 
CALL COMPARE (TEMP, 1, 8. CODE , IGOOD) 
IF(CODE.EQ.0)GOTO 1 


NINPTS = IGOOD 
NINPP1 = IGOOD + l 
ARRAN ERR ECHO NINPTS RAKKKEKKKAAKKK 


WRITE (*,2076)NINPTS 


CARRARKRKKKRERKKAKRE MODIFY NINPTS IF NEEDED KRARRRERRKKRKK 


16 WRITE(*, 2077) 
READ (*,2190)ANSWER 
e. TE (ANSWER. .EQ.'N!' .OR.ANSWER.EQ.! ee me ee 
IF ANSWER.EQ. '¥' .OR-ANSWER. EQ.'y')GOTO 15 


GOT 
A 17 CONTINUE : 
cee SKIP COST FUNCTION ENTRY IF NUMBER OF CONTROLS .GT. 1 RARKKKK 
IF(NINPTS .GT. 1) THEN 
GNSKED = 3 


wpTeoore 340 

ENDIF 
CRAEKKKAAAKKAEKKAKKRKERKRAKRRKEKRKERERERKRKRRREKRRERRRKRKRRRKRKKRREKKKRRREREREKRARKEKRRRKEKK 
CRARAARARRARRRARE ENTER THE NUMBER OF TIME INTERVALS = *xxxxkkRRKKR 


CRAEKKKRKAKKAKAKKKKKKRKAKKAKRKKAKRRERAKRKERERERERKERRKERRERERERERERERERERKRERERKERKKK 


20 WRITE (4 2080) 
(*,*)NSTAGE 
TF (NSTAGE .GT. 1000)GOTO 20 
NSTP1 = NSTAGE + 1 
EE (CHNGH EQ. 1970 780 
IF(FINAL .EQ. 1)GOTO 1520 


CRAKRKAKRKKKRRKKKRKARKKERKKKRKEKKKARKRKKRKRRKERKKKERKEKERRKKKKRKRKARRKKKRKKRRERKRRRKK 
CREKKAAKAKAKRKKAKREKRA INPUT THE MATRIX ARKKAKKKKAAKKKK 
CRAKKKAKRKRKKKKAKKKAKKARKKRKKERKKKERAKRKRERRKREREARKKRRKKRKERKERERRKKRRRRERRRKRKRERKRRRKK 


30 LOOP = 0 
WRITE (*, , 2090 
READ (*,2070)T 
CALL COMPARE (TEMP, ae on CODE , IGOOD) 
TFC CODE: steed sp Gete 
OPTION GOO 
GOTO (40,50 0) OPTION 


GOTO 3 
40 Wien c 100) 
50 WRITE(*,2110) 
GOTO 80 
60 WRITE(*, 2120) 
80 DO 90 I = 1,0RDERN 


DO 90 J = 1,ORDERN 
IF £0. J), THEN 
IF(OPTION .EQ. 1 9(I,J} = 1.0 
IF(OPTION .EQ. 3 x io) = 0.0 
IF(OPTION .EO. 2) THEN 
WRITE (#, 2130)1 J 
READ (*,*)Q(1I,J) 
ENDIF 
ELSE 
eae) = 0.0 
90 CONTINUE 
CRERAKAKRRAAKAKKA ECHO THE Q MATRIX KRAKKAKKKKRA 
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100 CONTINUE 
WRITE (*, 2140) 
DO 110 I=1,ORDERN 

110 WRITE(*, 2150) (9(1, J), J=1, ORDERN) 
IF(LOOP .EQ. 1) GOTO 30 


BARA AAAKAES MODIFY THE Q MATRIX IF NEEDED KAKAKKKKKRRR 


120 WRITE(*,2160) 
READ (*,2070)TEMP 
CALL COMPARE (TEMP,1,3,CODE,IGOOD) 
IF(CODE. 2E9.0) 0 )GOTO 120 
OPTION = 
GOTO(130, EO) 160) OPTION 
_ GOTO 120 _ 


a ee CHANGE ONE ELEMENT OF THE 0 MATRIX AKA RER 
130 se Zi 


Oo 
od 


I. .GT.ORDERN .OR. J.LT.1 .OR. J.GT.ORDERN)GOTO 130 
J) 
DERN 
Q(I,J),J=1,ORDERN) 
150 WRITE(* 2180 
READ (*,2190)ANSWER 


TF (ANSWER. .EQ.'N'.OR.ANSWER. EQ. 'n BET 160 
IF(ANSWER.EQ.'Y' .OR.ANSWER.EO.'y')GOTO 130 


150 
160 IF (FINAL.E0.1)GOTO Pozo 


CRAKKKARAKRKKAKRKAKAKRKAKRKRKAKAKKERAKRKERARAKRAKRARAKRRRKARARKRRKRKRARKRARARRRRRARKRA 
CRAKRKKAKKKAKKKKRKAKK INPUT THE H MATRIX RAAKAKKKKKKA 
CRARKAKRKAAKKAAKKAAKKAKAKKAKRKKAKRKKAKKRKRKAKRKKAKKRKARKRKAKRARAKRAKKAKKRAKRKRAAKKAKRKRKKERRKA 


D 1 
140 eee x” 2180} 


170 LOOP = 0 
Se a So 
READ (*,2070)TEMP 
CALL COMPARE (TEMP , i 3, CODE, IGOOD) 
~—IF(CODE. can aueore 70 
OPTION = IG00 
cet eo, 190 300) OPTION 


GOTO 1 

180 WRITE(*, 2210) 
c 210 

190 WRITE (* 2220) 
GOTO 210 

200 WRITE(*, 2230) 


210 DO 220 I = 1,ORDERN 
DO 220 J = 1,0RDERN 


220 CONTINUE 
CARAAKAKRARARRARA ECHO THE H MATRIX KAKAKKKKKRER 


230 CONTINUE 
WRITE(*, 2250) 
DO 240 I=1,ORDERN 


ey 


240 WRITE (* 2150) (H(T, 3). J=1 , ORDERN) 
IF(LOOP .EQ. 1) G 170 


C aki ac MODIFY THE H MATRIX IF NEEDED KRAKKKKKKKAKKK 


250 Gere eae 
D (*,2070)TEMP 
CALL COMPARE (TEMP , 1,3,CODE,IGOOD) 
IF(CODE.E ey geie 250 
OPTION = 
GOTO(260, 130 9390) OPTION 
GOTO 250 


Okan ARR CHANGE ONE ELEMENT OF THE H MATRIX KRKEKKKKKKKKK 
Cc 
260 WRITE saa 
= READ (*,%) 


oO 
~~ 


- 


f 


D 
270 eS te ,J) ,J=1,ORDERN ) 
280 WRITE(*, 2180 


READ (*,2190)ANSWER 
TE CANSHER, .EQ.'N! .OR. ANSWER. EQ. 2 290 
IF(ANSWER.EQ.'Y!'.OR.ANSWER.EO. GOTO 260 


GOTO 
290 IF(FINAL. “FO. 1)GOTO 1520 


C 

CRARAKKAARRRAKKRK RK KKK RRA ARKERAKK RRR AKA RRA AKA AK AK KERARKRRRER ARK RERARERRK RR 
CRAKRRRRKARKAK KKK INPUT R KAKRKKKRAKKK 
CRARKAKAKKARKKKRRAKKRARAKKKAAKKK RKEAKK RKRAKKKRARAKREKAERAKARARAKKKRRARKRRR 


300 WRITE(*,2260) 
READ (*,*)R 


c 
CRAKKAAKKAKKKAK KR ECHO R KAKKKKKKAKKA 


310 WRITE(*,2270)R 
CRAKKKKKKKAKKAKKKKK MODIFY R IF NEEDED KAEKKEKKKAKKEE 


320 WRITE(*, 2280 

READ (*,2190)ANSWER 
EE CANSHER. .EQ.'N' .OR.ANSWER.EQ. eae 330 
IF(ANSWER.EQ. 'Y!' .OR.ANSWER.EO.'y!')GOTO 300 


GOTO 
330 IF(FINAL. EO. "1 )GoTO 1520 
CRAKKKKKKKKKKKKAKRKKKKKAKAKKKKAKKAAKKAKE RK KAKKRKERERERRAKKRRERRRERKRARKK KK 


CAKKKKKRKKKKKKAKRKA CHOOSE TO ENTER EITHER A KRAKKKKKKAKKK 
CRHARKKRKKKKAKKKKKE CONTINUOUS TIME SYSTEM OR A KAEKKKKKKAKKKKEE 
CRKKKKRAKRKKKKKEAKKE TE TIME KREEKAKKKKKRKKK 


DISCRE SYSTEM 
CRKEKAKKKKKKKAKKKAKKKKAKKKKKKKAAAAKAKAKKAKAAAKKAAAKAKAKKAKKKKKKKKKAEKARKRREREKER 


340 WRITE (* 2230), 
(* 2070) TEMP 
CALL SOMe eee 0, 1, CODE, IGOOD) 
IF(CODE. -E9.0) .9)GOTO 3 40 
SYSTEM GOO 


TF (SYSTEM) 350 350, 360 
350 WRITE (*, 230 
READ (*, 2190) ANS WER 
IF (ANSWER.EQ.'Y'.OR.ANSWER.EQ.'y')GOTO 370 
SYSTEM = 1 
360 WRITE(*, 2310 
READ (%, 2190) AN WER 
IF (ANSWER.EOQ. ' TOR. ANSWER.EQ. 'y! )GOTO 590 
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SYSTEM = 0 
Geror 350 


e 

CHAKRA AKAKKAAKKRARKRKKKKARAKKAKRERKARKAKKAAKRKKARKARKKKAKARARKAKKKARKAKRKRERRRRRRAKRKA RAK 
CRAKKKRKKAKKAKKKKKKRK INPUT THE A MATRIX KRARAKKKRERKKEE 
CRAKKAKAKKKAKRKAKKAKKKKAKKAKAKKAKRKARERERERERERERKARRARRRRERERARRRKARERERERRERERERRE 


370 WR TEs, 2320) 
DO I=1 ,ORDERN 
BO 380 J=1 , ORDERN 
oO ee ,2330)1I,J 


R *)A(I,J) 
380 CONTINUE 


CAKKKKKKRKRKKKKKKKE DO NOT ALLOW CHANGES TO A and B KRAKKKKRRKRKRKKK 
See eAAKAAREX TF A DISCRETE TIME SYSTEM WAS ENTERED Ras NA a 


390 CONTINUE 
IF (SYSTEM FPS RE 


WRITE 305 
READ * 3070 TEMP 
GOTO 1520 
ENDIF 
CAKKKAKKKKKKKKKKKK ECHO THE A MATRIX KRAKKAKARARRRAAKR 
C 
WRITE (*, 2340) 
DO 400 I=1,ORDERN 
400 WRITE(*,2150) (A(I,J) ,J=1,ORDERN) 
BERRA AR ARERR MODIFY THE A MATRIX IF NEEDED KAAKRAARAKRARRR 
410 ee nee 
REA 2070 ) TEMP 


CATT. SOMBARE (TEMP, 1,3,CODE, IGOOD) 
IF (CODE. -£9.0) .0)GOTO 0 410 
OPTION = 
GOTO(420, 370 9450) OPTION 
GOTO 410 


CURE A Hn HH CHANGE ONE ELEMENT OF THE A MATRIX  *xkxkxkkRAAARAK 


420 WRITE (* Dsvaltalp 
READ (*,*)I,J 
TF(T PT 1 .OR. 1.GT.ORDERN .OR. J.LT.1 .OR. J.GT.ORDERN)GOTO 420 


D I= 
430 mes. 3180) TACT, J),J=1,ORDERN) 
440 WRITE(*. 2180 

READ (*,2190)ANSWER 
LE (ANSE ER.EQ.'N'.OR.ANSWER. EQ. 'n boos 450 
IF(ANSWER.EO.'Y! .OR.ANSWER.EQ.'y!')GOTO 420 


GOTO 440 
450 IF(FINAL.EQ.1)GOTO 480 


C . 

CRERKKRKKKEKRKKKRKRKRKRKKRKRKRKRKKREKRERAARKAKRRAKARAKAKARKAKRKAKKRKAKKKRKARRKRKRARKAKAARRRARERKAK 
CARERKRKKKKKKKKKEKKR INPUT THE B MATRIX KARR RRR 
CAKKKAKAKKKKKKAKAKARKAKKKAKAAAKKRRAKAAKKAKKRKRRARKRKRKARKKAKKRRRERKKRARKAAAKKARRRRARARRRKERK 


460 ee 2350) 
DO T=1 ,ORDERN 
DO. ‘470 ij = 1,NINPTS 
eae ,2360)I, J 
READ (*,*)B(I, iD 
470 CONTINUE 


— oo ECHO THE B MATRIX KAKRKKAKRAKKAKKK 
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480 CONTINUE 
WRITE (*, 2370) 
DO 490 I=1,ORDERN 
490 WRITE(*,2150)(B(1I,J) ,J=1,NINPTS ) 


Ce MODIFY THE B MATRIX IF NEEDED AARKKKAAKKKKKK 


500 WRITE (* 2160), 
READ (*, 2070) TEMP 
CALL aero cues 1,3,CODE, IGOOD) 
IF (CODE. -E9.0) .9) GOTO 500° 


OPTION = 
GOTO(510, 460, 9540) OPTION 
GOTO 500 
Ci sci ieki ii CHANGE ONE ELEMENT OF THE B MATRIX KAREKKAAKKREAKE 
G ig : ’ 
510 eT eae . 
READ (ae 2 
TFL ET 1 -OR. I.GT.ORDERN .OR. J.LT.1 .OR. J.GT.NINPTS)GOTO 510 
WRITE(*,2360)I,J 
READ (*,*)B(I,J) 
WRITE(*, 2370) 
DO 520 I=1,ORDERN 
520 ie 2150)(B(I,J) ,J=1,NINPTS) 
530 WRITE(*,2180 
READ (*,2190)ANSWER 


iF ANSHER -EQ.'N' .OR.ANSWER.EQ.'n oe 540 
IF (ANSWER.EQ.'Y!' .OR.ANSWER.EQ.'y!')GOTO 510 


540 IF (FINAL. "EO. 1)GOTO 1520 


e 

CARKKKKKKKKKAKRKAKKKRKKRKRKARERKRKKRRKRRERRERERKRRRRRRERRKRKRKRRRRERERERKRKRKRRKRKKRKKKKKKK 
CHAKKKKKKKKAKAKKKKK INPUT THE SAMPLE TIME....DT KRKEEKKKKRKKKKK 
CRAKKKAKKAKKKKAKKAKKKKKKKRRKKREKRKRRRKKRAKRKRKRRKKRKERERERKRERKERERKERRRKRRAKRKRKERKKAKK 


550 pees re 2380) 


2" )DT 
, DTFLAG = l 
CHAXKAKRAAAARRRAK TF A DISCHEDE) DEMEmo vac etinwoomen ce RED RAR ARARA AA 
CAAKARAAKAKKEAAKK = €=AND NO VALUE FOR Di HAS BEEN ENTERED TAs cae 
CRAKKKAAKK KA RKEA RR THEN PRINT OUT A MESSAGE KRRKAKRKAKRIK 


560 TEMRTREAG. ¥ 0385 QO) THEN 
READ (* SOTO VTEME 


GOTO 1520 
ENDIF 
CHAKKKKKKKKKKKAAE ECHO THE SAMPLE TIME....DT KRREKKEKKARKKAKE 
WRITE(*,2390)DT 
Ra ARR eee MODIFY THE SAMPLE TIME IF NEEDED KREEKKKAKKKAKKE 
570 WREATE (7 as 
READ (*,2190)ANSWER 
“TE (ANSWER. -EQ.'N'.OR.ANSWER.EQ.'n' ooo 580 
BREED EO.'Y' .OR.ANSWER.EO.'y')GOTO 550 


C 

CAKKRAKKKKKKKKRKEKKKAEKKEKKEKREKEKKRRERRRERKRREEKRERERERERERERREEKRRRRREREREKK 
CARKKAKAKKKKKKAKKKA CONVERT A and B TO) PHI and DEL RAKKKKKKKKK 
CoRR Ca A a ACI SI eR CEC De TCE CTC Te TC Fo ese ve cere 


580 IF(SYSTEM aero THEN 
CALL PHI EL(DT, ORDERN ,NINPTS ) 


ENDIF 
IF(FINAL.EQ.1)GOTO 1520 
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GOTO 780 


CARKRAKAKKAKAAKAKKKAKKKAKKKKKKKKKKKKKKKKKRRKRKRKRKRKRKKRKRKRRRRKRKRRRRRRRRRKRKRRRRRKRRRRKRKRKA 
CARAKKKKKAKKAKAKKKE INPUT THE PHI MATRIX KAEKKKKKAKKE 
CAKKAAAKAAKKKKARKKKKKKKKKKKKKKKRKRKRKKRKRRRRRRRKRKRKRRRRRRRRRRRRRRRRRRRRRRRRRRRKRRER 


590 CONTINUE 
DTFLAG = 0 
600 WRITE(*, 2410) 
DO 610 I=1,ORDERN 
DO 610 J= i, ORDE one 
eee ae 242 O)I 
READ (% *)PHI(I, "3) 
610 CONTINUE 


C 
CAARAAKARAAAAAAAK DO NOT ALLOW CHANGES TO PHI and DEL = *xxxAAAAAAAK 
CHAAAKKRAAAAAAAAA TF A CONTINUOUS TIME SYSTEM WAS ENTERED ****xxRAAAAK 


620 CONTINUE 
TF (S7s0EM EQ. OEY 
WRITE 
READ (%* * SOT STEP 
GOTO 1520 
ENDIF 


CRARAAKKRAAKKKRAA ECHO THE PHI MATRIX KAKKKKKKRKRK 


WRITE(*,2430) 
630 WRITE(*, 3150) (PHI (I, J) , J=1,ORDERN) 


C 
CRAKKAAKKKKAAAKKKR MODIFY THE PHI MATRIX IF NEEDED KREKKREAKKAKKER 


640 eee 2160) 
D (x , 2070) TEMP 
CALL COMPARE (TEMP, 1,3,CODE, IGOOD) 
IF (CODE. E SCOTS 640 
OPTION = G00 
GOTO(650, 600, 680) OPTION 
GOTO 640 


CRIA RRA RRA AAKA CHANGE ONE EISSN ON aA Ie IG UME SIG eee tele acta a aase 
650 WRITE(*, ero). 
REA (x, 


*)I,d | 
TF(I. LT.1 .OR. I.GT.ORDERN .OR. J.LT.1 .OR. J.GT.ORDERN)GOTO 650 


WRITE(* ,2420)I, J 

READ (*, /®)PHI (I, J) 

WRITE (*, , 2430) 

DO 660 I=1,ORDE 
660 eee ,2150 (PHT (I, J) ,J=1,ORDERN) 
670 WRITE 2180 

READ (*,2190)ANSWER 


TE CaNSHER. -EQ.'N'.OR.ANSWER.EQ. 'n ae 680 
IF(ANSWER.EQ.'Y! .OR.ANSWER.EQ.'y')GOTO 650 


GOTO 670 
680 IF(FINAL.EQ.1)GOTO 710 


CRAKKAAKRAKRAKKRAAKKRAKRAKKAKRAAKRAAKRAKRAAKAAKKRKARKARERAKRRARRK ARR AKRAARRKKR 
CRAAAKKAAKKRAARKRR INPUT THE DEL MATRIX KRAKKKKAAKKA 
CRAAKKRAAKRAKKAAKRKARKRAARKKARKARRARRAKRRARKKARRRARKAAKRARKRAAKRAAKRAKRAAKRAKARKK 


690 Na 2440) 
DO Ot ORDERN 


c 
WRITE (%, ,2450)1 


READ (* *)DEL(I, 5) 
700 CONTINUE 


CHRIGER A AHMAR AI ECHO THE DEL MATRIX KRRARKKKAKKRA 
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710 CONTINUE 
WALTENGA 2460 ) 
enn20 = 1 , ORDERN 
720 WRITE(*, 2150) (DEL(I, J), J=1,NINPTS) 


me a MODIFY THE DEL MATRIX IF NEEDED  **xxxxxAAAAR 


730 WRITE (4 2160), 
(* 2070) TEMP 
CALL ereeucecat 1,3,CODE, IGOOD) 
IF(CODE. 2£9.9) 0) GOTO 730 
OPTION = 
GOTO(740, 350 9770) OPTION 
GOTO 730 


Choco CHANGE ONE ELEMENT OF THE DEL MATRIX “xxkkkARAAAAA 


DO 750 I=1,ORDERN 
750 ae ,2150) (DEL(1,J) ,J=1,NINPTS) 
760 WRITE(*, 2180 
READ (*,2190)ANSWER 
Hy Eas -EQ.'N'.OR.ANSWER.EQ.'n! on 770 
IF (ANSWER.EO.'Y'.OR.ANSWER.EO.'y')GOTO 740 


GOTO 760 
770 IF(FINAL.EQ.1)GOTO 1520 


g 

CRAKKKKKKAAARKAKKKAKRRAERAAERAKARERAKRRERKAERRAERRERRERERRERRERRERERRARERERERRRKERKE 
CAARAAAAKARARAAAK WRITE ALL CURRENT INFORMATION TO THE pees aS eas st 
CAHARKARKKKKARKKKARKA OUTPUT FILE kaKAKKKKKKKKK 
CRARKARKARKARKKAKRKRARKERERAARRAERERARERERRAERARRERRARRERRRERKAERERRERKARRKRKRKA 


780 FINAL = 0 
WRITE (9,2030 
WRITE (9,2470 )ORDERN 
WRITE(9,2475 )NINPTS 

~—IF(GNSKED -EQ. 3) GOTO 805 
WaTTE (3 2480 NSTAGE 
WRITE (9, oe: 


.0 

DO790 I= 1 /ORDERN 

TRACEO = CEO +0 Q(I,1) 

790 WRITE (3 2490) (9 J), J=1, ORDERN) 
WRITE(9 228 250) 


DO DERN 
800 WRITE(9, 2450) GC, J) ,J=1,ORDERN) 
WRITE(9, 2270)R 
805 IF (SYSTEM) 810,810,840 
810 WRITE(9,2340) 
DO 820 I=1,ORDERN 
820 WRITE(9, zB ,J),J=1,ORDERN) 
WRITE (9, 2370 





DO 830 I=1,ORDERN 
830 Wee meee ,J),J=1,NINPTS) 
840 RE 2430) ° 
850 I= 1 ,ORDERN 


850 WRITE(S, ‘secon ,J),J=1,ORDERN) 
WRITE(9 , 246 
DO 860 I=1,0R 
860 WRITE (94 2480) (DEL(T, J) ,J=1,NINPTS) 
WRITE(9, 2030) 
IF(GNSKED .EQ. 3) THEN 
CRARARAKK NO OPTIMAL GAINS ARE TO BE CALCULATED KRARARARKRKK 
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GOTO 1010 
Dy 


N 
SE Eee 870,880 
870 WRITE(9,2500 
PNAME2= san dima TERMINAL STATES Control)' 


Ne 33 
GO TO 890 
880 WRITE (9, "9810) 
PNAME2= aut ns a over ALL STAGES)! 
PNAM2ZL= 30 
ieee A Ge RN HARARE RRR AKI Hick hhh heck kkk 
CRAARKAAAAKKKAA AKA INITIALIZE MATRICES PRIOR TO KAKAKKAKKKAK 
CAKKKKKKKKKKRAKKKKA ULATING OPTIMA KRAKKKKKKKAKKK 


CALC L GAINS 
CRAKKAKAAKKKKKKKKKKKKKKKKKKAAKKKKKAKKKKRKKKARKRKKKAKAKKKKKAKKKKKAKARKKRKRAKAKRK KKK 


a —s 


890 pore 
Bee T=1 , ORDERN 


oa = a0 
FM(I) = 0.0 
DO 900 J=1,ORDERN 
a a = 0.0 
“3 = 0.0 


B00 P(I,J)=h(0, J 
penne DO YOU WANT TO SEE THE GAINS TABLE ON THE SCREEN ? *****x*xx 


READ (*/ 2190) 
* 2190) ANSWER 


READ 
fae ER.EO. NOR ANSWER. 20.7 n- Beceen = 0 
: IF (ANSWER.EQ.'Y'.OR.ANSWER.EQ.'y')SCREEN = 1 
CAAKKKKKKKKKKKKKK PRINT HEADING FOR OUTPUT TABLE KAEKKKKKKEKER 
CAKKKKKAKKKKKKKKKK OPTIMAL GAINS KKEKEKKKKKKKKER 
IF(SCREEN .EQ. 1)TH 
WRITE (%, 2 20) (HDG), ,1=1,ORDERN) 
WRITE 2030 
NDIF 
ea ees 
WRITE(9,2030 
CAKKKKEKKAKKKKKKKKKKKRKKKERKKKKEKRKKKRKERKERKKRKERKKKRKAKKKRKKKKRKERKRKERKKKKRKEKKKK 
CARAKAAAKAKAKKAKK LOOP TO ITERATE THE RICATI EQUATIONS ARAKKKKAKKKAK 


CAKRKKKKKKKKKKKKRRRRRKRKAKKARKKRKARAKKKAKKAKAKKKKEKRKAKKKKKKKRKKKKKKKKKKKKKKKKKKKK 


DO 1000 KK=1,NSTAGE 
KREAL = NSTP1 =- KK 


DEN=0 .0 
DO 910 I=1,ORDERN 
DO 910 J=1,ORDERN 
910 EM(I) = EM(T) + + DEL (3, 1) * P(J,I) 
DO 930 ERN 
DO 920 oat ORDERN 
920 FM(I) = FM(T) + EM(J) * {BHI(J, I) 
930 DEN = DEN + EM (1) * DEL(I,1) 
DEN = DEN + 


E 
CxxkkRAAAA ENSURE THAT THE DENOMINATOR DOES NOT GO TO ZERO **xxRRAAAKR 


IF( DEN .E O )THEN 
WRITE (* 2S apes =e 
WRITE 9) 2530) KK=1 
NSTAGE = KK - 1 
GOTO 1007 

ENDIF 


c 
ce aera CALCULATE OPTIMAL GAINS FOR THIS STEP *****xxxXX 
DO 940 I=1,ORDERN 


143 


-FM(I)/DEN 
FTRAN(I) 


FTRAN(I) 
FNEG(KK,I) 
FM(I 
940 EM(I) = 0.0 
CARKARAKARARARKAR A PRINT OPTIMAL GAINS FOR THIS STEP RARKKAAKAKA 
IF(SCREEN .EQ. 1)THEN 
IF(ORDERN .GT. 4) THEN 
WRITE(*,2540)KK, KREAL, (FTRAN(I) , I=1,ORDERN) 
5 WRITE (*, 2541) KK, KREAL, (FTRAN(I), I=1, ORDERN) 
ENDIF | 
IF(ORDERN .GT. 4) THEN 
WRITE (9,2540)KK, KREAL, (FTRAN(I) , I=1,ORDERN) 


S 
SRETEN? ,2541)KK,KREAL, (FTRAN(I) , I=1,ORDERN) 


c 
CARKKAARARAKAAAKKRKA CALCULATE PSI(K,I,J) ARAEAAKRKKKAKRK 
e 
DO 950 I=1, oo 
DO 950 ‘J=1 ERN 
950 PSI(I,J) = PHT(L, Sy + DEL(I,1) * FTRAN(J) 
CARAKAKKKRKAAKAKKAKK CALCULATE P Cr lew) KRAAKKKAKKKKRE 
e 


DO 960 I=1,ORDERN 
DO 960 J=1,ORDERN 
DO 960 L=1,ORDERN 
960 GM(T, J) = GM(I,J) + PSI(L,I) * P(L,J) 
0 980 I=1,ORDERN 
” 50 980 J=1,ORDERN 
DO 970 L=1,ORDERN 
970 HM(I,J) = HM(I,J) + GM(I,L) * PSI(L,J) 
I,J.) = HM(I,J) + Q(I,J) + R * FTRAN(I) * FTRAN(J) 
980 HM(I,J) = 0.0 
DO 990 I=1,ORDERN 
DO 990 J=1,ORDERN 
990 GM(I,J) = 0.0 


CAA RAA RA DISCRETE TIME VECTOR FOR PLOTTING GAINS **x&xxkAARKX 
VTIME(KK) = 

“ 1000 CONTINUE 

CAAKAAKRAKKRARKA PO YOU WANT TO SEE THE GAINS PLOTTED ? AAKAKKKKARAK 

“1001 WR 2545) 


(* 2190) ANSWER 
"TE (ANSWER. .EQ.'N' .OR.ANSWER.EQ. 'n! }g0TO 1006 


IF(ANSWER.EQ.'Y'.OR.ANSWER.EQ.'y' )GOTO 1002 

‘a GOTO 1001 

CAKKKKAKKAKRKKARKKK LOOP TO PLOT OUT THE GAINS KRAKKKKKKAKKREAA 

G 

1002 DO 1005 GAIN = 1,ORDERN 

CRAKKKAKKKKAKKAKAKAKKK SET THE GAIN PLOT TITLE KRARAKAAKKKKAKKE 

c 
IF(GAIN.EQ.1)PNAME1L = 'FEEDBACK GAIN (Fl) FOR STATE X1' 
IF(GAIN.EQ.2)PNAME1 = 'FEEDBACK GAIN (F2) FOR STATE X2' 
IF(GAIN.EQ.3)PNAME1 = 'FEEDBACK GAIN (F3) FOR STATE X3' 
IF(GAIN.EQ.4)PNAME1 = 'FEEDBACK GAIN (F4) FOR STATE X4' 
IF(GAIN.EQ.5)PNAME1L = 'FEEDBACK GAIN (F5) FOR STATE X5' 
IF(GAIN.EQ.6)PNAME1 = 'FEEDBACK GAIN (F6) FOR STATE X6' 
Ie GAIN EO.8 PNAME1 = 'FEEDBACK GAIN (F7) FOR STATE X7' 
IF(GAIN.EQ.8)PNAME1 = 'FEEDBACK GAIN (F8) FOR STATE X8'! 
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PNAMIL = 31. 


c 
CRAKKKAKKKAKAKKKKAKRKKAAKKKKKAAKKEKKKKKKEAKKEKRKKKKKRKKKRARKKKEAKKKKKKKAKKKKAK 


CARKKKKKKKKKKKKKKR PLOT 88 GRAPHICS KAKRKKKKKKKK 
CARKAKRKKAKKKKKAKKKKKKKKKKKKKKKKKKKKRKKKRKKKKKKKKKKRRKRKKRERKERKKKK KKK 


C 
CRAAKRRARRAKKEAKKA SET sg INITIAL PARAMETERS FOR GAIN PLOT AAKKKKKKER A 
0 


BEGTIM = 
FINTIM = NSTAGE 
NPTS = NSTAGE 
DO Wsst 3) = a i 
1003 VTIMSS(J) = ((FINTIM - BEGTIM)/6.)*(J-1) 


ce 
eee GENERATE GAIN VECTOR FOR PLOTTING GAINS ***xAxxKKAAKK 
a DO 1004 KREAL = 1,NSTAGE ° = ' 
KK = NSTP1 - KREAL 
VY(KREAL) = FNEG(KK,GAIN) 
Boerne TEST LINE ee SELECTING PROPER COLUMN OF GAINS FOLLOWS oe 


CRAKKRAKKK E LL51 FOR COMPILED VERSION 

e WRITE(%, ae GAIN, KREAL, KK, FNEG(KK,GAIN) 

pe CONTINUE 

CARKAKKKAKAKKAKKKKK MAKE THE GAIN PLOT KREKKKAKAKK 
c 


IF (GAIN Hye 1) PAUSE 
CALL GRAPH( 99,991) 


C 
nnn he IS A HARDCOPY OF THE GAIN PLOT DESIRED ? **xxxxkKKAKK 


ea: 2595) 
READ (*,2190)ANSWER 
IF (ANSWER. EQ.'Y'.OR.ANSWER.EQ.'y') 


- CALL GRAPH(IOPORT ,MODEL, 1) 
CONTINUE 
1005 CONTINUE 
C 
CAARKKKKKAKAKAKKAAKAKKK DO YOU WANT TO CHANGE NSTAGE ? ARAKKKAKAKKAAKK 


g 
1006 CHNGN = 0 
WRITE (*, ea 
-~—READ (*,2190)ANSWER 
IF (ANSWER. EQ.'Y' .OR.ANSWER.EQ.'y' )THEN 
CHNGN = 1 


GOTO 20 
ELSEIF (ANSWER.EQ.'N'.OR.ANSWER.EQ.'n' )THEN 
GOTO 1007 
LSE 
GOTO1006 
DIF 


C 
CARKKKKKKKAKKKAKRK IS A PHASE PLANE DESIRED ? KREKKKKAKKRAER 


C 
1007 WRITE (* 2547) 
ee 2190) ANSWER 
eee CAN SWER.EQ.'Y'.OR.ANSWER.EQ.'y') THEN 
NSTP1 = NSTAGE + 1 
PHASE = 1 
GOTO 1025 


ENDIF 
IF (ANSWER.EQ.'N'!.OR.ANSWER.EQ.'n')GOTO 1010 
GOTO 1007 


C 
CAAKAKAKAKAKKKKKK IS A TIME RESPONSE DESIRED ? KARKRARKARKRKKKK 


e 
1010 PHASE = 0 

WRITE (* ,2550 

READ (*,2190)ANSWER 

IF(ANSWER.EQ.'Y!' .OR.ANSWER.EQ.'y!)GOTO 1020 


145 


PPAR LT RR eC CTIA LPI sg feCa 1510 
CRKKKKKKKKKRKAKKKKK GRAPH IS TO BE A TIME RESPONSE KEAKKKKKKKKKK 
1020 NSTP1 = NSTAGE + 1 
PLTYPE = 3 


c 
CAKKKKKKKKAKKKAKKK HOW MANY SECONDS ? AKAKKKAKKKKKKK 


C 
1025 WRITE (* 2560) 
(x *)TFINAL 


Cia INPUT DT IF NOT ALREADY KNOWN KAKKAKKKAAKK 
IF(DTFLAG . x Psaq) 0) THEN 
WRITE 


READ ba es po 


| eka CALCULATE FINAL VALUE OF K KEKEKAKKKKKKKE 
g 
KFINAL = NINA) 
TFTEMP = KFINAL 
IF(TFTEMP .LT. ne THEN 
KFINAL = KFINAL i il 
TFINAL = KFINAL * DT 
. ENDIF 
CAAKAKKAKKAAAKKK A ENSURE THAT ENOUGH GAINS ARE CALCULATED *xx*&xxxxxxx 
CRKKAKKKKKKKAKKAAK TO COVER THE DESIRED TIME RANGE KAEKKKKAKKKKKK 
€ 
IF (er Nae .EQ. 3) GOTO 1029 
IF((KFINAL- 1) ~GT. mS NCE) THEN 
MAXTIM DT * NSTAGE 
WRITE (* , 2561 )MAXTIM 
GOTO 025 
ENDIF 
CLERH KARE EI READ IN THE INITIAL STATE VECTOR KAAEKKKKKKAKK 


C 
1029 WRITE(*, 2565) 
DO 1030 I=1,ORDERN 
me ware (4, /2566)I 
READ (*,*)XKO(T,1) 
1030 CONTINUE 


rn one READ IN THE COMMAND INPUT VECTOR KAAKKKAKKKA 


WRITE (*, 2570) 
DO 1035 I=1,ORDERN 
ase ,2580)I1 
READ */*) TEUE (2d 
1035 CONTINUE 


C 
CAAKKKKKAAKK WRITE INITIAL STATE AND COMMAND INPUT VECTOR **xxxxxxxx 
CAAKKKAKKKK TO OUTPUT FILE KRKEAKKKAKKKK 
WRITE(9, 2568) 
WRITE(9, 2584 
DO 1036 I = 


WRITE(9, S585) T I, XKO(I,1),INPUT(I,1) 
1036 CONTINUE 


¢ 

CRAKKKKKKKKKKKKKK CHOOSE EITHER STEADY STATE GAINS (1 KKKKKKKKKKKK 
CHRKKAKAKKKAKKKKAKKR OR DYNAMIC GAINS 2 KREAKKKKKKKAKK 
CAHKKKKKKAKAKAKKAKKKK OR USER DEFINED GAINS 3 KEKKKKKKKAKK 
CAKKAKKKAKKKKKKAKK IF ONLY ONE CONTROL INPUT IS USED KKEEKKKKKKKKK 


c 
1040 ee Ome) ae 
IF (GAINCH NE. 3) THEN 
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ae ca 

READ (*,2070)TEMP 

CALL COMPARE (TEMP, 1,3,CODE,IGOOD) 
Th CODE. enn a RGEO 1040 

GNSKED GOO 


LSE 
GAINCH = 1 
DLE 


SE 
GNSKED = 3 
sara 
O(141 | i42, 143) GNSKED 
DER RRAAL RAR ANER USE STEADY STATE GAINS RARAKKAKAKKKKKE 
141 PNAME3 = Toor ice STEADY STATE GAIN SCHEDULE' 
PNAM3L = 34. 
.GOTO 1054 | = 
CRARKKEKKKKKAKKAKKK USE DYNAMIC GAINS KREKEARKKAKKKKKRK 
142 PNAME3 = 'OPTIMUM DYNAMIC GAIN SCHEDULE' 
PNAM3L = 29. 


GOTO 1054 
CAAAAKAKKARAARAAA IMPLEMENT USER DEFINED FEEDBACK GAINS **xxxxxAKAAK 
143 PNAME2 = ‘Implementing! 
PNAM2L = 12. 


1 
PNAME3 i a DEFINED GAINS' 


PNAM3L 1 
ie vanes .EQ. 1 ) GOTO 1043 
c F( 1 NAL NEO: 1 .AND. GNSKD3 .EQ. 1 ) GOTO 1043 
Cha RAKES AAKAKRKK INPUT USER DEFINED FEEDBACK GAINS RRRKKKRAKK KK 
1044 DO 1045 fr = 1,NINPTS 


B 2s), 
x oo) USERGN(L, J) 
1045 CONTINUE 


GNSKD3 = 1 
CARAKKKAKRARAKARK ECHO USER DEFINED FEEDBACK MATRIX RRRKKKAKKEE A 


C 
1043 WRITE(*, 2593) 
DO 1046 I=1,NINPTS 
1046 WRITE(*,2594)(USERGN(I,J), J=1,ORDERN) 


c 
CAAKKAKKARKKRKKRAA MODIFY THE USER DEFINED GAINS IF NEEDED ***xx*xxx 


C 
1047 WRITE (* 2160) 
(* 2070) TEMP 
CALL. COMPARE(TEMP ,1,3,CODE,IGOOD) 
IF(CODE. 2£9.0) .0)GOTO 1047 
OPTION GOO 
GOTO(1048, i944. 1052)OPTION 
GOTO 1047 


c 
CxXXAARRAR CHANGE ONE ELEMENT OF USER DEFINED GAIN MATRIX *x*kkxkARAKKK 


C 
1048 WRITE (4 2470) 
IF (I. LD. 1 Gensel .CIMNINPTS WOR. GeLT.1 
OR. J GT.ORDERN)GOTO 1048 
* WRITE * 2592)I 
READ (x « ADUSERANCT, a) 
WRITE 
D ech I=1,NINPTS 
1049 WRITE(?, 2594) (USERGN(I,J) ,J=1,ORDERN) 
1051 WRITE(*. 2180 
READ (*,2190)ANSWER 
TE (ANSWER. .EQ.'N'.OR.ANSWER.EQ. 'n! eos 1052 
IF (ANSWER.EO.'Y'.OR.ANSWER.EO.'y')GOTO 1048 
GOTO 1051 


é 
Cabin k ks WRITE USER DEFINED GAIN VECTOR TO OUTPUT FILE = *****xxxxx% 
1052 CONTINUE 
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WRITE( 9; 2595 
DO 1053 I = 1,NINPT 
WRITE(9, 2594) (USERGN(L, J) ,J=1,ORDERN) 
1053 CONTINUE 
WRITE(9, ek 
1054 aN EQ.1)GOTO 1520 


WRITES; 2593) 


IF -EQ. 0) GOTO 1050 
CRARKKRERAE ER ERE CALCULATE STATES FOR PHASE PLANE ARAKKKKAKKKAKRK 
CALL STCALC(1, Rene) 
CHARRKAKAKKAKR RA SET UP INITIAL PARAMETERS FOR THE RARKKKAKKKK 
ee PHASE PLANE PLOT RAKKAKKKKKKK 
NPTS = KFINAL 
FORK He PLOT THE PHASE PLANE KAKKKKKAKKAKK 
C = 
- ce CALL GRAPH(99,99,2) 
CRAAKAKKAKAAKRAKKR IS A HARDCOPY OF THE PLOT DESIRED 2?  *xxxxkRAAKKK 
C 


Ae see 
READ (*,2190)ANSWER 


IF (ANSWER. EQ.'Y' .OR.ANSWER.EQ.'y!)CALL GRAPH(IOPORT,MODEL, 2) 


CONTINUE 
‘ GOTO 1010 
CRAAAAKARARKRAKKA DO YOU WANT TO SEE THE TIME RESPONSE RAKKKRRKAKRKAK 
CRAKKKKKKRARKKKRA TABLE ON THE SCREEN ? KAKKKKKKKKKA 


C 
1050 WRITE(*,2591) 
READ (*,2190)ANSWER 
TF (ANSWER BQ. (Ne Seah Qual SCREEN 0 
IF(ANSWER.EQ.'Y'.OR.ANSWER.EO.'y!)SCREEN = 1 


C 
CRAAKAAAKRARARKAK SELECT HOW THE STATES ARE TO BE PLOTTED ***xAxAKKAAKK 


151 WRITE(*, 2598 
READ (*,2070)TEMP 
CALL COMPARE(TEMP, L, 13) CODE , IGOOD) 
IF(CODE. :EQ.0)GOTO 1 


STPLOT 
C 
CAAKKARAKKAKAKKKKA LOOP TO PLOT OUT STATE TRAJECTORIES RRKRKKKKKARRK 
C 
a DO 1500 STVAR = 1,ORDERN 
CRAKKAKKKAKKKKRARKKR IS THIS STATE TO BE PLOTTED ? RAKRKRKRKRRRRKRRE 
C 
IF(STPLOT .EO. 2) T 
WRITE (*, 2 as STAR . 
READ (*,2190) ANSWER 
ee  EOv IiaOR . See e ee Rn eienas = 0 
IF (ANSWER.E .'Y! (OR. ANSWER.EO. 'y! PLOTCH = 1 
ENDIF 
Curie fevededsictoiedeieledc terest ici fete i fed H HRM RII SH EHO PERO EH TF 
CREKKKKRKRKRKEAKKKKRKKE PRINT: HEADING FOR OUTPUT TABLE KRRERKRRKRRKEREE 


CRAKKKAKKKKKKKKKAKK TIME RESPONSE KRRAEKRKRARRKAKKRE 
CRAKKKRRKRKKEKRKRKRRKRKKEKERRRRRERERRRKRRRRRRERRRERRRRRERRRRRRKRERRRKRRERRRRERERRRERK 


IF(STVAR .E yy 


TH 
IF (SCREEN 2901 
WRITE(*,2 25) (uae (2), , 1=1,ORDERN) 


WRITE (*, 2030 
ENDIF 
WRITE (9, ,2525) (HDG2(I) , I=1, ORDERN) 
WRITE(9, 2030 
; DIF | 
CRAAKAKKK SKIP PLOTTING IF NO PLOT IS DESIRED RAAKKKKKKRK 
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CxARKAKKR = BUT MUST CALCULATE STATES ON FIRST TIME THROUGH ****xxxxxxx 


g 
IF(STVAR .NE.1 ) THEN 
TF (STPL OT .EQ. 2 GOTO1499 
IF(STPLOT .EQ. .AND. PLOTCH .EQ. 0) GOTO1499 
, ENDIF 
CAAKKKAKKKKKAKKKKKK SET THE PLOT TITLE BASED ON THE : KEEKKKKKKKKEE 
CAKKKKKKKKKKKKKKK STATE SELECTED KEAEKKKKKRKKKK 
e 
IF(STVAR.EQ.1)PNAME1 = 'X1l TIME RESPONSE' 
IF STVAREO.3 PNAME1] = 'X2 TIME RESPONSE' 
IF(STVAR.EQ.3)PNAME1] = 'X3 TIME RESPONSE' 
IF(STVAR.EQ.4)PNAME1 = 'X4 TIME RESPONSE’ 
IF(STVAR.EQ.5)PNAME1 = 'X5 TIME RESPONSE’ 
. IF(STVAR.EQ.6)PNAME1 = 'X6 TIME RESPONSE' = 
- IF(STVAR.EQ.7)PNAME1 = 'X7 TIME RESPONSE' 
IF(STVAR.EQ.8)PNAME1 = 'X8 TIME RESPONSE’! 
PNAMIL = 16. 


C 
CHAKKKKKKKAKKKAKKKKKKKKKRKEKKERERKEKKKKKKEKKEKREKEREKRKRKERKERKEKKRRERERRRRKRRRERERE 


CRAKAKAKKAAKAKAAX CALL SUBROUTINE TO CALCULATE THE STATES ****xxxxkkxX 
CRA KKAKKKKKKAKKAKKKKKAKK RAK KKAKKKA KKK KKAKRAKKAKKKKAKRERKRRAKRARRRARERRRRKR 


c 
CALL STCALC(0,XKO,STVAR, SCREEN ) 


C 
CRARAAAAAAAAARARK = SKIP PLOTTING IF NO PLOT IS DESIRED = *#*AAAAAAAKK 


C 
PE eeeE oe .EQ. 3) GOTO1499 
IF(STPLOT .EQ. -AND. PLOTCH .EQ. 0) GOTO1499 é 
CRAKKKKKKKKKKRKAKKKKKKAKKKRKRKERKKRKRKKKKEKAKKKKRKRKRREKRRKREREARRKRRKRRRRKRRRREK 
CAKKKKKKKKKKKKRKRRE PLOT 88 GRAPHICS KRAEKKAKKKKKK 
CHKKKKKKRAKKKKKKKKKRKKKRKEKRKRRAKKRKRKEKKKRKRKRKKKRKRKRERKRKERKERRAKRRKRKRRKKKRRKERERKRRK 
Gia ieikihk SET UP INITIAL PARAMETERS FOR THE RAAKKKAKKKR 
CARKKKKKKKKKAKKKK STATE TRAJECTORY PLOT KREKKKKKKKK 
1055 BEGTIM = 0.0 
FINTIM = TFINAL 
< NPTS © = KFINAL 
DO 1060 J = 7 
VYSS(J = INPUT(STVAR,1) 
-— VTIMSS(J) = ((FINTIM - BEGTIM)/6.)*(J-1) 
1060 CONTINUE 

RT eax ck a PLOT THE STATE TRAJECTORY KEAEKKKKKKKK 
C 

IF(STVAR . aeeNG ) PAUSE 
. CALL GRAPH(99,99,3) 
CAKKAAKKKKAKAKKKKK IS A HARDCOPY OF THE PLOT DESIRED ? KEEKKKKKKKKK 
c 

Be x” 37303 

READ (*,2190)ANSWER 

IF (ANSWER. EQ.'Y'.OR.ANSWER.EQ.'y' )CALL GRAPH( IOPORT, MODEL,3) 
CONTINUE 


1499 CONTINUE 
1500 CONTINUE 


e 
CAAKAAAAAAAAAKA PRINT OUT THE AVERAGE VALUES OF ALL STATES *****xxxxxx 


C 

WRITE(*,2030 

WRITE(X, P2596 

WRITE (9, 2596 

DO 1505 I=1,ORDERN 
WRITE (* meer ae ee eo 
cere 9,2597) I,AVG(I),1I,AVG2(1I),1I,MAXVAL(I 

1505 CONTINUE 
Res 2030) 
PAUSE 
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~» 


C 
CRARKAKAAAAK AAR IS ANOTHER RUN OF OPTCON DESIRED ? AAKKKAKKARKKR 
C 


1510 WRITE( 


2600) 

READ (*,2190)ANSWER 

ze (ANSWER -EQ< |: Gages BOs 2 aa 1520 
IF (ANSWER.EQ.'N'.OR.ANSWER.EQ.'n')GOTO 1530 
GOTO 1510 


WRITE (9, 2030) 
(x 


CRRA ae PRINT MENU OF OPTIONS AAKKKKKKKAKK 


C 
a2 0 


C 
CRAKARRARKAEAA 


1530 


WRITE(*, 2610) 
READ (*, 2070)TEMP 
CALL COMPARE(TEMP,1,11,CODE,IGOOD) 
TF (CODE .EQ .0)GOTO 1520 
-OPTION = IGOOD 
IF(OPTION .LE. 4) THEN 
IF(NINPTS .GT. 1)THEN 
WRITE(*, 2620) 
READ(*,2070) TEMP 
GOTO 1520 
ENDIF 
ENDIF 
IF(OPTION .EQ.2 .OR. OPTION .EQ.3) LOOP = 1 
IF(OPTION .E0.8) GAINCH =1 
IF (OPTION, -EQ.10 .AND. GAINCH .EQ.1) GAINCH =2 
GOTO(20, 230,100, 310,390,560,620,1040,10,780,1510)OPTION 
GOTO 1520 


STOP 


CRERKAKKRRARK 


C 


G 
CAKKKKAAKAKKAAKAKKKKAKKKKAAKKKKAKAKKKKAKKKRKKKKAKKKRKAARKKRKKKRKARAKKKAKKKAKKKAAKKKRR 


CRARKKAKKAAAKARAKKR FORMAT STATEMENTS KAAKKKKKAKKKKK 
CARAKKKAKKKAAKKKKAKAKKKAKKKKAAKAAKKAAKKAAKKKKRREAKKARAKKKRARAKKAKKKAAKKKAAKKKAKAKAKRK : 


G 
2000 


2010 


Zuo 


FORMAT(/,5X,'OPTCON minimizes the following cost', 
+! :' //,5X,'JIJ = MIN ( X''(N) * H * X(N) + ', 

+! Sum/( ee * 9 A X(k) + Ul (k)!, | 

+' * R * U(K)))',//,5X,'The output of the program is the', 

+' feedback eae matrix, F transpose (FY) a / Sx, (which, when', 
+' multiplie by the State Vector (x eve K,'y1lelds a scalar’, 
+' control.(U).',///,5%,'The following recursive equations ", 
+'were derived using dynamic pro eam Ag «eee 

Onna aean at the terminal time.(N) and working backwards:',//) 
+'(1) Fil (k) a PERRI 1 DBI a ae + R)',3X, 


'FI1(Q)=0! ; ; 
Bute ae = PHI + DEL*F''(k)!,27X 'PST(0)=0! ,/ 8X 
+'(3) P(k) | = PSI! '(k)*P(k-1)*PSI(k) + O + F'' (KS *RAF(k)', 4x, 
+! PCO) SE) a, 


FORMAT(/,5X,'You may enter a system with either single or', 

+' multiple control staraisi)) /7c eee ced Seen with only one', 

+' control signal is entered,',/,9x,' then the optimal gains can', 

+' be generated as described',/,9x,' above. These gains may then', 

+' be implemented into the',/,9x,' state equations to obtain a', 
time response of the system.',/,9x If you choose to enter a', 

+' system with multiple control signals, | ,/,9x,' then yee must!, 

+' enter the feedback gains manual i The user defined! ,/,9X, 

+' gains option exists BOn che Single control input system also.', 

+, /1// AX, First enter Ene wp LooLen . 

+1X,' identification ( NOT to exceed 20 characters ).',//, 

+10XK,.' PROBLEM “five... 3 

FORMAT (A20) 

FORMAT(! | 7 70 (=a) ee 

FORMAT(///,5X,'OPTIMAL CONTROL PROGRAM! ,/) 
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2050 
2055 


2100 
2110 
2120 
2140 


2150 
2160 


FORMAT (Se ,/,' PROBLEM IDENTIFICATION:' ,5X,A20, 

FORMAT(5X,//,' Select the type of printer that you are', 

+ using "Sve 
+! ( Answer ere 2 ) 'e//, 
ees EPSON or THINKJET! , / 
One. Z LASERJET! ,//, 
+tLOX, “ANSWER. GI. oc se ae 
FORMAT(5X,//,' Enter the ORDER of the system (up to 8). ', ) 
FORMAT (A2 
FORMAT (5X,//,' tw the NUMBER OF CONTROL INPUTS ae tO ce) 7s 
+ OX, | NOTE...NO OPTIMAL GAINS will be generated if you enter',/, 
+ ak more than | one control input',//, 
~ Bx  MANISWER.....5«* 
Four); 5x, tne NUMBER oF CONTROL INPUTS = ',I1) 
FORMAT (/,5X,' | areas changes to ayaa OF CONTROL INPUTS ? : 

swer or 
FORMAT (5X, if RAtee the’ NUMBER of TIME INTERVALS WN) ‘over', 
' which the cost function',/,' is to be! 
' minimized. (MUST NOT exceed 1000) 

FORMAT (10X, //, Doesmihe cost function (J) include the State’, 
+! TRAJECTORY over all Stages ?',/, 
+! ( Answer ie Or 3 ea! //+ 
is ye set . equal to pie IDENTToY Matric / 
+1lORee 2 YES. .Each diagonal element of Q will be entered! 
+! A ad - 
+1 OKpS NOs. . «SECC. 20 at to the ZERO Matrix .',//, 
TG ee NSWER a0 0 0 0 0 6 0 the 
FORMAT (9X, //,' The states are weighted equally for the', 
+' TRAJECTORY over all staqes. 
FORMAT(9X,//,' Enter the elements of the Q matrix.',/, 
+' (State ‘weighting matrix for TRAJECTORY over all stages)', /) 
FORMAT (9X,//,' The state TRAJECTORY is not included in your', 
+' cost function. 
FORMAT (6X, Boa ce ritarje=s |. ) 
FORMAT(//,5%, | The f aca ey) 
FORMAT (2X, 8(F8.3 rn 
FORMAT(//,5X, 'Do you want to change an element of the matrix?', 
+//,10X '1) "'YES...a SINGLE element. ' 
+#10X,'Z YES...the ENTIRE Matrix.'! fe 
+10X,'3 NO’ //, 
OA, PUN WER sce! cones 0 cereus 
FORMAT (/,5X,' Which element of the Matrix do you want to!', 

~+! Change. 2", 
+5X is the ROW and J is the COLUMN, .enter Lig —— 
FORMAT(10X,//,5X,' Any other changes? (Answer Vacreny ', ) 
FORMAT (Al) 
FORMAT(10X,//,' Does the Coot function $9) include TERMINAL', 
+' States ? ( Answer r 3 
t10X, 11) YES...Set H Saal” eomeHe sbenrrry Matrix .',/ 
+10X,'2) YES...Each diagonal element of H will be entered! 
+' separately ey . 
LOK 2) eNO a...«. SQ fags to the ZERO Matrix .',//, 
+10X, NWANSWER. +e ssseeeeee! 
FORMAT (9X, J//,' aa states, are weighted equally for the', 
+! TERMINAL states. 
FORMAT (9X,//, BRnee the elements of the H matrix. Poy es 
+' (State weighting matrix for TERMINAL states)' J) 
FORMAT (9X,//, The TERMINAL states are not included in your', 
oc Ose function. 
FORMAT(6X,'H(',I1, Pale) oa 
FORMAT(//,5X,' The ab Matrix ',/$ 
FORMAT //,5X,' Enter the value of the cae ee 
+5X, Control input SS anand EAeOr) ao) ass 2) 
FORMAT /,5X,' The scalar R 8.4 
FORMAT(/,5X,' Any changes to R" 5 (Answer yorn) ', ) 
cee 1/1, 4% 


+! ou want to read in the A and B matrices for a CONTINUOUS', 
+! E sysuen | aa 
BMH c. oi 'e) sl ate) o)c.6 oc ateminOne tenon Nese eeastenelenes 6 se cee es Enter Q', 
ay 4X, 1 Tf you want to enter. the “PHI ee DEL matrices for a', 
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2540 
2541 
2545 
2546 
2547 
2550 


2560 
2561 


2505 
2566 
2570 


2960 


fo DISCRETE TIME systen, ' 
+! 


»/ 48, 


+10X, (ANSWER... ss eeces. 


.-Bhiter o'>//, 


FORMAT(/, 5xX,' You will enter the A and B matrices. Se 
+S -[s this correct? 
FORMAT (/, 5x, | “You will enter the PHI and DEL matrices. ',/, 
+5X, | tae coo this correct ? 
FORMAT 5x, on ST the $j meee of the plant matrix...A.',/) 
FORMAT 6x, 'A(',I LL 
ee re No NaN are to A or B will be allowed PSGaus cmt’ 
or yee have entered a DISCRETE TIME sygcenie 
ENTER tO Conbinde........ cu 
FORMAT ip a: The A Matrix (Plant Matrix)', /) 
Praha Ske = the Bey) of the contro distribution', 
ma rix... 
FORMAT 6X, 'B , Le =! ) 
FORMAT(//, 5x : one B Neer (Control Distribution jiatrix)', h) 
FORMAT (5X, sK, “EAbee the SAMPLE INTERVAL.....DT = 
ee 5X,/,' No changes to DT will be allowed Bakaucet 
ae fae entered a DISCRETE TIME system! fy 
ENTER to Continues... ... 
FORTY ex q > The SAMPLE INTERVAL DT = '! non 4) 
FORMAT (/ ag! Any changes to the SAMPLE INTERVAL ? (Answer', 
FORM r( Sx, LL! 1 shtes the Reeements of the PHI matrix.',/) 
FORMAT (6X, 'PH Lee 
Ee SK / No Renee té PHI or DEL will be allowed because',/, 
5A ne have entered a CONTINUOUS TIME system) hs 
on fit ENTER to continiie. .ce..22.. 5 ee 
FORMAT //,5X%,' The PHI Matrix' ,/) 
FORMAT(5X,//,' Enter the elements of the DEL matrix.',/) 
FORMAT 6X, 'DEL(' 3 ve ) 
FORMAT(//,5X,' The bet aes Bay 
FORMAT 1//,5%, ' The ORDER of the system = '! ,T1) 
FORMAT(///,5X,' The NUMBER OF CONTROL INPUTS = ',T1) 
FORMAT(///,5X,' The NUMBER of TIME INTERVALS = ',13,/) 
FORMAT 2X,8(F8. 31) 
FORMAT(//,' Minimum TERMINAL STATES Control' ) 
FORMAT(//,' Minimization over ALL STAGES' ) 
FORMAT(// ,4X, ou want to see the gain schedule table on', 
+' the screen i gill 5X, aero V One) ae) 
aC ' NEG', AL! 
TIME! ,! RT IME 118, 4(A5,5X),/, 
STEP! ,, INDEX! ,T18,4(A5,5X),//) 
FORMAT(//," REAL' ,/, 
TIME REAL! ,T20, Algae oe ie 
' INDEX TIMES 'T20 4(A4,8X //) 
FORMAT(/ , Optimum gains. are reached after We 3: stages.' 
go The Padus is be nmenacee early in order t cm 
revent a division by z eel 
FO i, eae ae PL Lon V4 (EB. 4, 2X /,T16,4(F8.4,2X)) 
FORMAT(' '/2(14,2X) ,T16,4(F8.4, 2X) ) 
re I, 4X,’ Do you want to see the gains plotted ?', 
+//,5X ' (Answer Or Toa 
co Ce 4X,' Do you want to change the NUMBER OF STAGES ?', 
+//,5% (Answer ¥ Omi ee) 
FORMAT (// Pe MOE ae to see a PHASE PLANE of Xl .vs. 
+! (ee / 8x," Answer Or 1 mae 
FORMAT (//, ax, 7 er want to see a time response of your', 
+! seem ? PT ORR '(Answer yooron) ', 


FO mee ,4%,' For how many seconds ? ', 

Pe Sean 5x,//' The o timal gains are computed for only ',F8.4 , 
' secon Please enter a smaller number. ! 

FORMAT (5X, i peer. the -_ of the INITIAL STATE vector ', 


FORMAT 64, X'S 0 a) 
FORMAT (5X,/' enter the “elements of the COMMAND INPUT vector-R.',/) 


FORMAT(6X,'R(', 
FORMAT(TS,' N + 114, intrray STATE! ,T35,'COMMAND INPUT' ) 


FORMAT (Iv 1y i loe FO. 4, T37,F9.4 
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2590 


2529 
2600 


FORMAT(10X,//,' Select a gain schedule...( Answer 1,2,or 3 )',//, 
tLOX, -L) | UsewSTEADY STAT oe Al gains over all steps Bh 
TlOx, "2 Use DYNAMIC fate oll 

Eee STEADY STATE ee DEFINED gains .',//, 


FORMAT (//, 4X,' Do you want to see the a response table on', 
+' the screen ?'! 1) 3X X, ' (Answer ee or Ee 


FORMAT (6X,/,' eed GAIN F ( =? ', )j 

FORMAT //, 5X, ' The USER ae GAIN Meer ic af) 

FORMAT (' i 8(E7. 4,1X) 

ee il, 4X, " De you want a hardcop Of this paiot 7? *, 

nswer 

pemantc? ee Ok; AVERAGE AND MAX VALUES OF ALL STATES',//) 

FORMAT (6X, ' moerece Value ote eel ,E12. 4,/, 
6X,! Average Value of Tl, a8 = BE 2uiaey 
6X,' Maximum Value of X', ous = " B1d 4, 1/5 

FORMAT(//- 5X, 'Do you want to Pilots. oe ! 5 

+//,10% ig) ‘ALL state trajectories.',/ 


t1OR, 12) Only SELECTED state trajectories. my 
NO state Peaeccor ses - hs 
tILOX ANSWER. 5 00s 
a Sie! Do you want to see a PLOT for state el! 2 ey 
X,' (Answer orn 
FORMAT (//, aX, ' This concludes he optimal control program', 
+! (OPTCON) . oe ene you want to run the program', 


+' again? niger Vaca) oa 
2610 FORMAT(///,5X,' SELECT ONE OF THE FOLLOWING OPTIONS: bees 
TO 1 Biancememomlmmermr Of SLAGHS 2. ....200cccccceens No. 
+10X,/,! 2) Change the TERMINAL state weighting matrix..... A’, 
+10X,/,' 3) Change the TRAJECTORY state weighting eee rr Vo, 
re), 4) Changestme CONTROL weighting factor.....cceeeeeR',/, 
+10X,/,' 5) Change the present A and B matrices',/, 
Piven oO) Change themomirer INTERVAL. . 0.255550 cb essecees pt /, 
+10X,/,' 7) Change the present PHI and DEL matrices',/, 
+10X,/,' 8) Change (or select) different FEEDBACK GAINS | ae 
#10X,/,' 93) Input an ene NEW SYSTEM', 
+10X,/,' 10} NO CHANGES. 
+10X,/,' 11) EXIT the program’ oT, 
+10X, 'SELECTION...( MUST be a number Bebe lecncelly jas. .+. |, 
2620 FORMAT(5x,/,.' No. change EO ems Pacanece! ~s allowed ecause', 
' you have entered a MIMO system.' 
+ Sx Hit ENTER to continue......e.eee'!, ) 
END 
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APPENDIX C 
OPTCON SUBROUTINE LISTINGS 


The following subroutines are written in MICROSOFT Fortran and are to be 
used on an IBM compatible system. These subroutines are required. by the main 
OPTCON program found in Appendix B and by the PLOTS88 subroutine found in 
Appendix D. A brief synopsis of the subroutine functions is given below. 


= - 


PHIDEL - Convert the continuous time A and_B system matrices 
to the corresponding discrete time ® and I matrices. 


PROD - Perform simple matrix multiplication of two matrices. 
Maximum dimension of the matrices is limited to eight. 


SUM - Perform simple matrix addition or subtraction of two matrices. 
Maximum dimension of the matrices 1s limited to eight. 


COMPARE - Test a user input response to determine if the response 
lies within the range of allowable integers. 

CLRSCR - A DOS command that allows the monitor screen to be cleared 
prior to the generation of a new graph. 

GOTOXY - A DOS command that positions the cursor to a designated 


coordinate position on the monitor screen. 


STCALC - Calculates the time response of system by iterating the 
discrete state equations. 


SNOdebug 
E - LL63sub 1L2JULY87 
C OK SDL 
G NEW output format 
E for states 
CRERKKKKKAKKKKKKKKEKEEKRERKEEREREREREEERRRKKEKARRRERKEKEERRRKREERERRKRERERRER 
CRAARKKKRKKKRKRKKKRKKRE SUBROUTINES KRERKEKREKEKKRRA 
CREKAKKKKKEEKREREKREEKEEEKEKRERERKEEREREREEEERRKKEKREREREREREERERREREKRRRRKERRER 
C 
C 
. SUBROUTINE PHIDEL(T,ORDERN,M) 

COMMON /BLK1/ A,B,PHI,DEL 

INTEGER*2 ORDERN,I,J,ERFLAG 

REAL*4 A(8,8 Peet ey Meee 

+ PSIT(8,8) , TERM(8,8) ,NEXTRM(8,8) ,ARATIO(8,8), 

; + TRATIO,ERROR,K 

ERROR = 1.E-7 

ERFLAG = 
é TRATIO = T*T/2. 

DO 1 I 1 ,ORDERN 

O1J=1,0R 


J = 1,ORDERN 
TERM(I,J)= A(I,J) * TRATIO 
IF(I .EQ. THEN 
PSIT(I,1) = T + TERM(I,I) 
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LSE 
<a = TERM(I,J) 
1 CONTINUE 


‘ore) 


K = 2. 


2K=K +1. 
TRATIO = T/K 
DO 3 =, an 
pes J = ie 
ARAT O(I, oy AL, J) * TRATIO 
3 CONTINUE 


_CALL PROD(TERM,ARATIO,ORDERN, ORDERN, ORDERN ,NEXTRM) 


DO 4 I = 1,0ORDERN | 
DO 4 J = 1,ORDERN 
TF (ABS (NEXTRM(T_ J) ) eee ERROR) THEN 
Bae = ERFLAG ih 


ENDIi 
TERM (I,J) = NEXTRM(I,J) 
4 CONTINUE 


- 


— 


enn acl O) DEN 
DO 5 I = 1,ORDERN 
DO ee J = 1,ORDERN 
PSIT( I» i= THRMGiede)) + PSIT(I , J) 
5 CONTINUE 


ERFLAG = 0 

GOTO 2 
ENDIF 

CHKKKKARKKKKKAKKKKE NOTE THE DUAL USE OF ‘TERM! HERE KREKEAKKKKAKKAAKK 
CALL PROD(A,PSIT,ORDERN, ORDERN, ORDERN, TERM) 
DO 6 I = 1,ORDERN 
DO 6 J=1, OR THE 
EE a HEN 
ae = 1.+ TERM(I,I) 


SE 
= PHI(I,J) = TERM(I,J) 


NDIF 
6 CONTINUE 
: CALL PROD(PSIT,B,ORDERN, ORDERN,M, DEL) 
RETURN 
END 
C 
& 
GIGS SISSIES IOI ICSI ITSO ITI TODS IIIA IASI ID RASA A aaa bk 
e 
; SUBROUTINE SUM(M1,M2,0PER,N,M,MSUM) 
INTEGER*2 N,M,OPER,I 
: REAL*4 M1(8, Se 4368, 8) ,MSUM(8,8) 
DOUt on % ie 


=1, 
1 MSUM(T, °s= 0. 0 


CAAKRKKAAKAKAKKAAKKKK PQ YOU WANT TO ADD OR SUBTRACT 2 *XRAARKKKKAKKKKK 


F(OPER) 2,2,3 
ChARKERRERE RES xk SUBTRACT KRERKEARKRERKERKRERE 
2 DO 20 I =1,N 
DO 20 =e = 1,M . 
unr, J) = M1(I,J) - M2(I,J) 


i> 


20 CONTINUE 


GOTO 40 
CRAKKKAKKKKKAKKKKKK ADD REKKKKKKKKKKKK 


3 pO 30 IT=1,N 
DO 30 J=1.M 
MSUM(I,J) = M1(I,J) + M2(I,J) 
30 CONTINUE 
40 RETURN 
END 
C 
C 
Cdededeiedeiei died Aric tee ede vee a Oe Aede fei de ice cide fk Pci ACA 3 fo 
C 
SUBROUTINE PROD (M1,M2,ORDERN,M,L,MPROD) 
-INTEGER*2 ORDERN,M,L,1,J,K ie 
REAL*4 M1(8,8), 1308" '8),MPROD(8,8) 
DO 1 I=1,ORDERN 
DO 1 J=1,L 
1 MPROD(I,J)=0.0 
DO 2 I=1,ORDERN 
2 J=1,L 
pO 2K=1,M 
2 MPROD(I,J) = MPROD(I,J) + M1(I,K) * M2(K,J) 
RETURN 
END 
C 
C 
Saetedcdedeiedekicd dietetic ici iok AI RIEI AR ARIA A RII RR ARR 
C 
- SUBROUTINE COMPARE (TEMP, VALMIN, VALMAX, CODE, IGOOD) 


INTEGER*2 IGOOD , CODE , VALMAX , VALMIN 
CHARACTER*2 TEMP 


IGOOD = -1 
IF (TEMP .EQ.'0')IGOOD=0 
IF(TEMP.EQ.'1' )IGOOD=1 
IF(TEMP.EQ.'2')IGOOD=2 
IF(TEMP.EQ.'3')IGOOD=3 
ea: IF(TEMP. EQ. 14! ‘TGOOD=4 
IF(TEMP.EQ.'5' )IGOOD=5 
IF (TEMP .EQ.'6' )IGOOD=6 
IF(TEMP. EQ. 17! ‘STGOOD=7 
IF(TEMP.EQ. '8' ) IGOOD=8 
IF(TEMP.EQ.'9! )ITGOOD=9 
IF TENE EQ: |10, (ee 


a IF (TEMP .EQ.‘'11')IGOOD=11 
IF(IGOOD.EQ.-1 .OR. IGOOD.GT.VALMAX .OR. IGOOD.LT.VALMIN) THEN 
CODE = 9 
EL 
CODE = 1 
ENDIF 
RETURN 
END 
C 
c 
Stet i AAA Rca Tei Pie TA INET TE TIC FICCI CRAIC 
C 
‘ SUBROUTINE CLRSCR 


INTEGER*2 ne 
COULVALENCE (ot 26( C2,IC(2))9( C3, 1C(3)), (C4, 1¢Ca> 
DATA Ic/16 15) (C2. 16939 sehak, 
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C *** Write ees Code Ee at KKK 


MRbeeee. 1) C1,C2,C 
1 FORMAT (1X, 4A1) 

RETURN 

END 
c 
Cc 
Be ee ee Hea ARERR EE SHA AREEN RRR RRR ARAKI 
c 
P SUBROUTINE GOTOXY (ROW, COLUMN) 
: xARKK Position Cursor by Row,Column ****x* 

INTEGER*2 IC(4) ,ROW, COLUMN, L ; _ 

“ CHARACTER*1 C1,C2,C5,C8,LC(5) 

CHARACTER*S CBUFF 

EQUIVALENCE et SN cage 8a SEB TCKS)), 

+ CBUFF,LC(1 

DATA IC/16#1B,16#5B , 16#3B, 16#66/ 
4 L=10000+100*ROW+COLUMN 
C *** Write Esca : Codes to a Character Buffer *** 

WRITE(CB ae L 
- 2 FORMAT(IS) 
C *** Write Escape Codes to Display *** 

WRITE (*, 3 3c chee tc(2), 1¢(3) ,C5,LC(4) ,LC(S) ,C8 

3 FORMAT ( : 

SETuEN 

END 
S 
C 
ea eek A AAA AR RAGA ATA E REE AREA KAR Rios 
G 
, SUBROUTINE STCALC(PHASE, XKO, STVAR, SCREEN) 


i ae CALCULATE THE STATES ITERATIVELY *** 
RAK _X(k+1) = PHI * X(k) + DEL * U(k) *** 


COMMON /BLK1/ A,B,PHI,DEL 
COMMON /BLK3/ VIIME, VTIMSS, Wi vio, VAnKoS Vass 
COMMON /BLK4/ KFINAL, NSTAGE, NSTP1, ORDERN, , GNSKED , USERGN, FNEG, 


+ INPUT , DT, AVG, AVG2 ,MAXVAL ,NINPTS 
INTEGER*2 KFINAL ,NSTAGE ,NSTP1,ORDERN , GNSKED , STVAR, SCREEN, 
+ PHASE .M,OPER ,NINPTS ,NINPP1 
REAL*4 A(8,8) B(8, 2 2 , PHI (8, € By DEL (8. 2) ,VTIME(1002), 
+ VTIMSS (9) AM Viss (9). NEG G(1000,8), INPUT (8, oe 
+ Dr: PHIE cia ty. (BRLEQ(B 8) (D pees 8), ,ROWF(8,8), 
+ XKO(8,1),XK(8,1),XKP1 1(8 22, TxXSS 9) VXY¥SS(9), 
A DELING(8 2) ,AVG(8), 'AVG2(8) , STSUM, STSUM2 , MAXVAL(8), 
7 USERGN(8 
CAR ARERARAEA RE-INITIALIZE THE STATE VECTOR KRAEKKKKKKEKKK 
C 
© 5 J = 1,ORDERN 
XK(J,1) = XKO(J,1) 
5 CONTINUE 
RENE ARE ARAARAE IE RE-INITIALIZE THE AVERAGING SUMS KKRAEKKRKKKKAKK 
CAKKKAKAKKAAKKKAKAKKK AND MAXIMUM VALUE STATE VECTOR KAREKAKKKKKKE 
G 


STSUM 0.0 
STSUM2 0.0 
MAXVAL(STVAR) 0.0 


C ; 
CAAAKARARAAAKAKARK LOOP TO ITERATIVELY CALCULATE THE STATES ***xxxxxKAKX 
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DO 70 K = 1,KFINAL 
KPRIME 


= NSTP1 - K 
TIME = K * DT 
IF (PHASE .EQ. 1) THEN 
VY¥(K = XK(2,1 
VTIME(K) = XK(1,1 


XK (STVAR,1) 
VTIME (K) TIME 
CRARARKARARKAKKKK SUM FOR COMPUTING AVERAGE STATE VALUES RRAKKKKARARK 
STSUM = STSUM + XK(STVAR,1 
STSUM2 = STSUM2 + (XK(STVAR,1) * XK(STVAR,1)) 
CRAKKRAKERKARERKER SEARCH FOR MAXIMUM VALUE OF THE STATE **x&x&xxkRXKK 
IF( ABS( XK(STVAR,1) ) .GT. ABS{ MAXVAL(STVAR) )) 


+ MAXVAL(STVAR) = XK(STVAR,1) me 
Che ENDIF 
C 
CAAAAKAKKKKKKKAKK LOOP TO SELECT THE PROPER FEEDBACK GAIN *X*X*XxXx*xxkRKKK 
CRAKKKAKKKKKKKKKK ELEMENTS FOR THIS TIME STEP KEERKEKRKRKKRRE 
C 
DO 40 J = 1,ORDERN 
GOTO(10,20,35)GNSKED 
CRAKKKKKKKKKAKKKKE USE STEADY STATE GAINS (GNSKED=1 ) FXXXARAARRKK 
10 ROWE = FNEG(NSTAGE,J) 
OTO 
CRARKKKKKAKKKKKKKK USE DYNAMIC GAINS (GNSKED=2 ) XX XAAAAARAAK 
20 ' IF(K .LE. NSTAGE) THEN 
ROWF(1,J) = FNEG(KPRIME,J) 
ROWF(1,J) = 0.0 
ENDIF 
30 CONTINUE 
CAKKAKKKKKKKKKKEK USER DEFINED GAINS (GNSKED=3 ) XA&XAARAAARKK 
oo CONTINUE 
A 40 CONTINUE 
CRAKKKKKKKKAKKKKE PAD THE DEL AND ROWF MATRICES KAEEKAKEKAERE 
CRAKKKKKKKAKRKKKKE WITH ZEROS KRARKKKAKKE 
Ca IN ORDER TO MULTIPLY PROPERLY IN PROD AAKKKKKKKKK 
NINPP1 = NINPTS + 1 
ee DO 50 I = 1,0ORDERN 
DO 50 J = NINPP1,ORDERN 
“DELG a = 0.0 
ROWF(J,I) = 0.0 
: 50 CONTINUE 
CAARKKKAKKKAKKRKRKRAKRKRKREKRRERRARRERRERRRRKRRERRRKRERRRARRRRRRRRRRKRRRRRRRRKKRERKRE 
CRAKKKKKAKKKKKKKKK CALCULATE THE NEXT STATE K(k+1) KEEKKKARKRKRE 
CARKKAAKKAKRKRKKKKKKRKRKRKRKRKRKRKRRRKRKRERRKRKRKRRRRRRRERKRRKRKRRRRRERKRRKRKRKKRRRKRRRRRRRK 
Cc 
M=1 
IF(GNSKED .NE. 3) THEN 
CRAKAKKKAKKAKKKRKK USING OPTIMAL GAIN SCHEDULE KEEKKKKKKEES 
eee PROD( DEL, ROWF , ORDERN , ORDERN , ORDERN , DELROW) 
CRAKKKKAKKKKKKAKKK USING USER DEFINED GAIN MATRIX KREKRKEREKEEE 
CALL PROD(DEL,USERGN, ORDERN , ORDERN , ORDERN, DELROW) 
ENDIF 
OPER = 1 


CALL SUM(PHI,DELROW,OPER ,ORDERN ,ORDERN, PHIEQ) 
CALL Nee Cate Bie GHOSE ERC R EINE Bat 
CALL ae DELROW, INPUT, ORDERN, ORDERN ,M, DELINP) 


OPER 
CALL SUM( PHIE X,DELINP ,OPER,ORDERN ,M,XKP1 
CRAKKKKAKKKAKKRKKRKERKKAEKRKEREREKRKRERERRKKE RRR RERERRERERRRRAKRERERRAKRKRKKK 


e 
CARAAAKAKAAKARRRKA NEXT 29 LENES ARE TEST DINES: TO VEREFY "<x 5077% 
CHAKRAKKAKKKRARKARER PROPER CALCULATION OF THE STATES RERKKRKKK KKK 
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CAAKKAKAAKKKKAKAKKK FOR A SECOND ORDER OPTIMAL EXAMPLE KRAARKAKKKRARRAE 


1041 


ANNAANANANAANANAANANAAANN 


C2640 


PE iat oore K 
WRITE (*, 2615 
DO 1041 I = 1,ORDERN 

WRITE (*, 2620)DEL(I, 1),ROWF(1,1I),(DELROW(I,J),J=1,ORDERN) 
CONTINUE 


pelea 3 = ERN 

WRITE (%, 2630) (PHI(L, J), J=1,ORDERN) , ( 
J=1 /ORDERN), (PHIEO(I, J 
CONTINUE 


WRITE(*, 2635) 
DO 1043 I = 1,ORDERN 
WRITE(*, 3640) (PHIEO(I, J) ,J=1,ORDERN) , KK(T, De PROM? a 
_. CONTINUE . 
WRITE(*, 2645) 
DO 1044 I = 1,ORDERN 
WRITE(*, 3680) PHIEOX(I, 1) ,DELINP(I,1),XKP1(I,1) 
CONTINUE 
FORMAT (/, TIME STEP = ',I3,/) 
FORMAT(/, ae DEE eee cor TRAN', 744, 'DELROW' ) 
FORMAT 15, F10.4,T20,F10. 4, 735, BRLO ; 
FORMAT(/,T10,' PHI'.T38,' DELROW | ‘Tee 'PHIEQ ') 
FORMAT Peace 4), 8X, 9(F10, 4) 8X, 2(FL0 4)) 
' 51, 'PHIEQX') 


y gai ROW 


| 


FORUAI(/, 110; ae eae 
FORMAT 2(F10. 4) ,8X,F10.4 ex F10. M4) 


C2645 FORMAT he Ti03! "sue Re 10x, DEINE ', 12k," xXKP1, ') 

C2650 FORMAT 3(Fl yy 

eee ee ee tee a AA K AAANARAKAAARER AKER RRR RRR IR 
GRR An An Re ae NEAT #24 LINES ARE TEST LINES TO VERIFY KAKKKKKK 
CAAAARKARAKAKA PROPER CALCULATION OF THE STATES RAKRKKKK 


EAAARAAKAAAK FOR A FOURTH ORDER USER DEFINED GAINS EXAMPLE ****xxxx 
CRARAAAKKKAKRAAKKRKRAKAKK RR RARKKKKARRAAKK RK RK AAAKRKRRARRAKKKRRARAKKKKRRR 


e ; 

G WRITE eee K 

C WRITE(9,2615 

E DO 1041 I=1, 

C WRITE(9, 5620) (FHI (I. J),J=1,ORDERN), (DEL(I,J) ,J=1,NINPTS ) 
C1041 CONTINUE 

C WRITE(Y, 2625) 

C DO 1042 I 1,ORDERN 

g a WRITE(9, , 2030) (DEEROW(T. apa Sere aeee eta 

C USERGN(J,I),J=1,NINPTS 

C1042 CONTINUE 

C WRITE(S, 2635) 

E DO 1043 I 1,ORDERN 

C WRITE(9, ae ii oe J=1 CRDeRAT PHIEQX(I,1), 

& XKP 1 ( { XK(I pe 

C1043 CONTINUE 

C2614 FORMAT(/,' TIME aed a, 

C265 FORmT(/ 210,' PHI',T57, 'DEL' 

C2620 FORMAT(TS,4(F7.4,2X),T50, O(F7.4 nx) ) 

CZezommenmal(/,.110,' DELROW' ,T52, 'DELINP ' T67, 'USERGN' ) 

C2630 FORMAT(TS,4(F7.4,2X),T50,F7.4,T62, 2(F7.4, )) 

C2635 FORMAT(/,T10,' PHIEQ',TS1, ae x! ,T63, 'XKP1' ,T74, 'XK') 

C2640 FORMAT(T5,4(F7.4, 2X 750, E7.4 PALO, F7.4.T70,F7.4) 

coaisacncnaenstaek SEO PRINT ue ae STATE TABLE KREEKKKKKKRKKEA 
CRAKKKKKAKKKKKKAKAKK Y ONCE KEEKAKRRRAKKRE 
c 


IF(PHASE .NE. 1) THEN 
IF(STVAR .EO. 1) THEN 
IF(SCREEN .EQ. 1.) THEN 
IF(ORDERN .GT. 4 ) THEN 
WRITE (*,2670)K, TIME, (XK(I,1),I=1,ORDERN) 


E 
EE cect TIME (XK(I,1),I=1,ORDERN) 
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IF(ORDERN .GT. 4 ) THEN 
WRITE(9,2670)K, TIME, (XK(I,1),1=1,ORDERN) 


S 
fe ae ,1),I=1,ORDERN) 


END 
2670 Eee ',14,T7,F8.4,T15, a ee a jas, 4( E1024) oe 
2671 FORMAT | 14,77,F8.4,T15,4(F10.4,2X)) 
ENDIF 
ENDIF 
Re RRNA GET READY FOR THE NEXT ITERATION KEEKKKKAKEKE 
C ‘ 
DO 60 IT=1, ORDERN 
XK { pe ARP LC Le) 
60 CONTINUE 
- 70 CONTINUE es ; 
CRARKKAKAKKRKKRRKK CALCULATE THE AVERAGE OF THE STATE RAKKKKKKKKKK 
CRAEKKRKKKKAKKAKKKKKE BEING CONSIDERED ON THIS CALL KREKKKKRKAKRKKE 
e 


tora -NE. 0) THEN 
AVG(S TVAR) = STSUM/KFINAL 
AVG2(STVAR) = STSUM2/KFINAL 
ENDIF 
RETURN 
END 
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APPENDIX D 
PLOT88 GRAPHICS SUBROUTINE LISTING 


The following code is written in MICROSOFT Fortran and is intended to be 
used on an IBM compatible system. This graphics subroutine must be linked with the 
two program segments found in Appendices B and C. In addition, the Fortran, Math 
and PLOTS88 libraries must be linked. 


SNOdebug 
C LL63GR OK SDL 
S IZ etIE  o/ 
CRAKKKKKKKKKRKKKKKKKRARKRKERKARRKRKRERERRRRKARRRRRRRRRRERARRRRRRRRRRKRRRRKEKRRK 
CRAERKKKKKRKKKKKRKKK SUBROUTINES KRERKKRKKRKKRKR 
CRAEKRKKKKKRRKKKKERKRKKRKRRRKKRERERERRERRRKRERRRRERRRRERRRARERRERRRRRRRKRRRRRRRRR 
Cc . 
C 
: SUBROUTINE GRAPH(IOPORT ,MODEL,PLTYPE) 
IMPLICIT REAL*4 (A-Z) 
en /BLK2/ BEGTIM,FINTIM,NPTS, 
SNAML, YNAML , PNAMIL, "PNAM2L, PNAM3L 
” comion /BLK3/ VTIME, VTIMSS, VX, VYSS, VXXSS, VRYSS 
COMMON /BLKS/ XNAME,YNAME, PNAME1, PNAME2, PNAME3 
ae NPTS, IOPORT, MODEL, XNAML , YNAML, 
NCHAR1, NCHAR2, NCHAR3, PLTYPE ne 
_REALM4 VTIME(1002), WE(1002) ,SASL, YAXL ,VTIMSS(9) , VYSS(9), 
XORGN , YORGN , VXXSS(9 VXYSS(9), 
XLO, XHI, MEO rrL INCRMT 
” CHARACTER*30 XNAME , YNAME 
_CHARACTER*51 PNAME1, PNAME2 
IF (MODEL -EQ. 99) THEN 
CRRKRAKKK RAK RRR A SEND TO MONITOR KRERERKERRRKRKK 
CALL CLRSCR 
XORGN = 1. 
YORGN = 0.80 
Bea A AK EAR Een SEND TO PLOTTER KRRAERKREREREAERKR 
XORGN = 3.20 
YORGN = 1.76 
ENDIF 
10 CALL GOTOXY(10,25) 
WRITE (*,*) — 7 Plotting Data' 
(PLTYPE .EQ. 1) T 
Cha RARE RARE ERERER See THE GAINS KRREKRRKRKRKRRRERRE 
AAAL = 5. 
SOFF = 0.25 
XNAME = 'DISCRETE REAL TIME INDEX (k)' 
SNAML = - 
YNAME = 'GAIN TRAJECTORY' 
YNAML = 15 
PNAME3 = ' ! 
PNAM3L = 1 
ELSEIF (PLTYPE Ose 2 a 
CAAKKKKKAKKKKKKKKKR eee cae PHASE PLANE RAEERRREKRRRRRRRE 
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XAXL = 4.0 
XOFF = 0.29 
XNAME = 'X1 STATE' 
XNAML = - 
YNAME = 'X2 STATE' 
YNAML = 8 
PNAME1 = 'Xl vs. X2 PHASE PLANE' 
PNAMIL = 21 
XORGN = XORGN + 0.65 
CERRER AAR OOO PLOTTING THE TIME RESPONSE RKAKKKKKKKKKKK 
XAXL = 5.0 
XOFF = 0.25 
XNAME = 'REAL TIME (sec)' 
XNAML = -15 
YNAME = 'STATE TRAJECTORY ' 
ae YNAML = 16 a 
2 ENDIF 
YAXL = 4.0 
ASPRAT = 0.70 
CHARHT = 0.23 
CHRHT2 = 0.8 * CHARHT 
CHAKKKAKKKKKKKKKRE PLOT TITLE LOCATIONS REAEKKKKKKAKRKKKKE 
cea = oor + (XAXL-PNAMLL*ASPRAT*CHARHT) /2. 
Bre = ee + (XAXL-PNAM2L*ASPRAT*CHRHT2)/2. 
ane = XOFF + (XAXL-PNAM3L*ASPRAT*CHRHT2)/2. 
NCHAR1 = ifix(PNAMLL 
NCHAR2 = 1f£1x(PNAM2L 
- NCHAR3 = 1f£1x(PNAM3L 
CALL PLOTS(0,IOPORT,MODEL) 
CALL Dee en eee 
- CALL ASPECT(ASPRAT) 
CALL SCALE(VY, ae NPTS ,1) 
IF (PLTYPE EON THEN . 
C This scalin spplies when the X axis represents DISCRETE TIME 
CALL SCAL :* IME , XAXL, Roe 1) 
CALL STAXIS 15, .20, .12, .080,0 
- ELSEIF (PLTYPE .EQ. 2) T THEN 
E ee Season ap lies when the X axis represents a STATE 
XHI = VITME 1 
YLO = Viel 
YHI = VY(1 
DO 15 J = 2,NPTS 
F ( VTIME(J) .GT. XHI ) XHI = VTIME(J 
IF ( VTIME(J) .LT. XLO ) XLO = VTIME(J 
IF VY(J) .GT. YHI ) YHI = VY (J 
IF VY(J) .LT. YLO ) YLO = VY(J 
15 CONTINUE 


XRANGE = XHI - XLO 

YRANGE = YHI - YLO 

IF( YRANGE .LT. XRANGE ) THEN 
INCRMT 


= XRANGE/XAXL 
VY (NPTS+1) = YLO - ((YAXL*INCRMT - YRANGE)/2.) 
VTIME(NPTS+1) = XLO - INCRMT/2. 
INCRMT = XRANGE/ (XAXL-1. ) 
ELSE 
INCRMT = YRANGE/YAXL 
VTIME(NPTS+1) = XLO - ((XAXL*INCRMT - XRANGE)/2.) 
VY (NPTS+1) = YLO - INCRMT/2. 
INCRMT = YRANGE/(YAXL-1.) 
ENDIF 
VY (NPTS+2) = INCRMT 
VTIME (NPTS+2) = INCRMT 
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20 


oe Seo 1556 20,.12,.080, 2) 


This scalin 
VTIME 
VT IME 


hae +1 
NPTS+2 


= BEGT 


plies rag the X axis represents REAL TIME 
y= (VEIME (NPTS) ye eee 


CAE Sbaais(.l5,.20,.12, 


ENDIF 


BERS LA 
DELTAX 
LASTX 

Foun. L 
DELTAY 


ee ene 
VIMEGiPist2 

FIRSTX + DELTAX*XAXL 
UNE rot! 

Vp t2 


~—CALL SYMBOL 


LASTY FIRSTY + DELTAY*YAXL 

MP (PLIYPE .EO. 1 .OR. PLTIYPE .EQ. 3) THEN 
ey BEGIIM 
eae . (FINTIM - BEGTIM)/XAXL 


DO 20 J = 1,7 


(((LastE - FIRSTX)/6.) * (J-1) ) + FIRSTX 
(((LASTY - FIRSTY)/6.) * (J-1) ) + FIRSTY 


FIRSTX 
DELTAX 
BER Sas 
DELTAX 
Puno? y 
DELTAY 


= FIRSTY 

= DELTAY 

CALL PLOT Ne 
L 


CALL PLOT eae 5 
CALL PLOT(0.00,YAXL,2 

CALL AXIS(0.0,0.0,XNAME, XNAML,XAXL,0O.,FIRSTX,DELTAX) 
CALL STAXIS(. SZ oF wea So, 2 

CALL AXIS(0.,0. AME, YNAML , YAXL 790. ,FIRSTY , DELTAY) 
CALL Sinsor era? el ks, CHARHT , PNAME1, 0. ,NCHAR1 


<i 

a 

oS 

Wn 

Ww -~ 
rn, 
C1 cy 
eee Cy 


CONTINUE 
Vi riso 
VIIMSS 
VAXS 
VXASS 
Vass 


Ww 
WwW OW Oo-~™ . 


CALL) SYMBOLGETAZ,PTYZ,CHRHTZ ,PNAMEZ ,0. ,NCHAR2 
PIAS; PIYS , GnenaZ, , PNAME3 ,O.,NCHAR3 


CALL LINE(VTIME,VY,NPTS,1,0,0) 
PEGeen Se 0 ) THEN 
IF( LASTY.GE.0 )CALL CURVE(VTIMSS,VYSS,7,-0.1) 


ENDIF 
IF(PLTYPE .E ee ae 
IF ( FIRSTX On ane 
IF ( LASTS. GE.0 yCALL CURVE (VXXSS,VXYSS,7,-0.1) 


ENDIF 
ENDIF 
Garr PLOT (O..,0.,999) 


RETURN 
END 
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